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CHAPTER  1 


INTRODUCTION 


1.2  General  Comments 


Adhesive  bonding  Is  increasingly  used  to  fasten  (metallic  to 
metallic  or  metallic  to  composite)  structural  components  together. 
This  is  because  in  many  present  day  applications,  conventional 
fasteners  such  as  bolts,  rivets,  welds  etc.,  are  unsuitable, 
especially  if  the  components  are  made  of  polymeric  or  composite 
materials.  Penetration  methods  (i.e.,  drilling  holes,  etc.)  cause 
high  stress  concentrations  and,  in  the  case  of  composites,  sever  the 
fiber  reinforcement  which  in  turn  reduces  the  strength  of  the 
joint.  On  the  other  hand,  bonded  joints  tend  to  be  damage-tolerant 
due  to  the  high  damping  behavior  of  the  adhesive  layer  and  less 
expensive  due  to  lower  fabrication  cost.  The  use  of  adhesives 
increases  the  joint  strength,  distributes  the  loads  more  evenly,  and 
enables  alternative  jointing  methods  to  be  reduced  or  eliminated. 
Dissimilar  materials  (e.g.,  steel,  aluminum,  plastics,  glass,  etc.) 
can  be  joined  together  by  bonding  even  where  it  is  impossible  to  gain 
access  to  either  side  of  the  joint,  thereby  increasing  the  design 
flexibility. 

Adhesives  are  not  free  of  disadvantages,  however.  Adhesives  are 
polymers  and  as  such  have  time  dependent  (viscoelastic)  moduli  and 
strength  properties  which  are  susceptible  to  environmental  effects, 
especially  temperature  and  moisture.  Most  polymeric  adhesives  are 
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rate  sensitive  materials  and  hence  exhibit  viscoelasticity. 
Furthermore,  certain  types  of  epoxy  resins  have  been  found  to  be 
nonlinearly  viscoelastic  in  character.  The  nonlinear  viscoelastic 
behavior  is  typified  by  a  stress-enhanced  creep.  Basically,  at 
elevated  stresses  the  material  moduli  seem  to  soften  and  the  creep 
progresses  at  accelerated  rates.  Time  dependent  properties  of 
adhesives  raise  serious  questions  regarding  their  long  term 
reliability  or  durability  under  creep  or  other  more  complicated 
loading  conditions.  A  delayed  failure  (creep  rupture)  long  after  the 
initial  design  and  fabrication  process  is  possible.  Thus,  methods 
are  needed  by  which  long  term  failures  on  the  order  of  a  structure 
design  life  time  (perhaps  as  long  as  5-20  years)  can  be  anticipated 
and  thereby  avoided.  Such  a  process  is  especially  important  in 
applications  where  failures  may  be  life  threatening  as  is  the  case 
for  automobiles,  airplanes,  missiles,  etc. 

Failure  in  an  adhesive  joint  can  occur  in  one  of  two  ways:  (i) 
adhesive  failures  that  occur  at  the  interfaces  between  the  adhesive 
and  adherends,  and  (ii)  cohesive  failures,  which  occur  either  in  the 
adhesive  or  in  the  adherends.  The  determination  of  the  strength, 
failure  and  reliability  of  an  adhesive  joint  requires  both  an 
understanding  of  the  mechanisms  of  adhesion  and  a  knowledge  of 
deformation  and  stresses  in  the  joint.  The  mechanisms  of  adhesion  is 
closely  related  to  the  chemical  and  physical  properties  of  the 
adhesive  polymers.  The  deformation  and  stress  states  can  be 


determined  once  the  geometry,  loading,  boundary  conditions  and 
mechanical  properties  of  the  constituent  materials  of  the  joint  are 
known.  The  mechanical  properties  of  the  adhesive  and  adherend 
materials  enter  the  stress  analysis  via  constitutive  models,  which 
relate  strains,  temperature  and  moisture  gradients  and  density  to 
stresses  and  fluxes  in  the  joint.  The  chemical,  physical  and 
mechanics  aspects  of  the  constituent  materials  enable  the  formulation 
of  appropriate  constitutive  models  for  adhesive  joints.  The 
determination  of  stresses  allow  the  prediction  of  the  strength, 
failure  and  reliability.  In  macromechanics  sense,  of  adhesive  joints. 

The  stresses  in  an  adhesive  joint  depend  on  the  geometry, 
boundary  conditions,  the  mechanical  properties  of  the  regions 
Involved,  and  the  type  and  distribution  of  loads  acting  on  the 
joint.  In  practice,  most  adhesives  exhibit,  depending  on  the  stress 
levels,  nonlinear-viscoelastic  behavior,  and  the  adherends  exhibit 
elasto-plastic  behavior.  Most  theoretical  studies  conducted  to  date 
on  the  stress  analysis  of  adhesively  bonded  joints  have  made 
simplifying  assumptions  of  linear  and  elastic  and/or  viscoelastic 
behavior  in  the  interest  of  obtaining  closed  form  solutions. 

A  good  understanding  of  the  process  of  adhesion  from  the 
mechanics  view  point  and  the  predictive  capability  for  structural 
failures  associated  with  adhesive  bonding  require  realistic 
theoretical  analysis  methods  to  determine  stress  distributions  in  the 
joint.  The  finite  element  method  is  the  most  powerful  analysis  tool 
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that  can  be  used  to  determine  stress  and  displacement  fields  in 
complicated  structures. 

At  present  there  are  numerous  computer  programs  available  for 
analyzing  bonded  joints.  However,  most  of  these  computer  codes 
incorporate  linearly  elastic  material  behavior,  and  some  allow  for 
nonlinearly  elastic  and  plastic  behavior.  Computer  programs  which 
incorporate  viscoelastic  material  behavior  are  quite  often  limited  to 
the  simple  spring-dashpot  type  of  model  for  linear  materials.  Such 
inaccurate  modelling  of  the  constitutive  behavior  of  the  structure 
can  seriously  compromise  the  accuracy  of  the  analytical  predictions. 


1.2  Objectives  of  Present  Research 

The  primary  objective  of  the  present  research  is  to  present  a 
finite-element  computational  procedure  for  the  accurate  analysis  of 
adhesively  bonded  joints.  With  this  aim  in  mind,  a  nonlinear 
viscoelastic  analysis  code  (NOVA  for  short)  has  been  developed.  The 
finite  element  program  NOVA  uses  linearly  elastic  elements  to  model 
the  adherends.  The  adherends  may  be  represented  as  isotropic, 
orthotropic  or  laminated  composite  materials.  The  large 
displacements  and  rotations  experienced  by  the  adherends  in  many 
types  of  loading  are  accounted  for  by  the  updated  Lagrangian 
description  of  motion  presented  in  Section  2.1.  It  should  be  noted 
that  this  description  is  valid  only  for  small  strains. 

The  adhesive  layer  is  modeled  using  a  special  element  that 
employs  a  multi-axial  extension  of  Schapery's  nonlinear  single 
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Integral  stress-strain  law  as  the  constitutive  equation.  The  finite 
element  formulation  for  the  viscoelastic  material  representation  has 
been  described  in  detail  in  Chapter  3.  The  element  library  contains 
an  eight  noded  isoparametric  element  which  employs  quadratic 
interpolation  functions  to  represent  the  displacement  field  as  well 
as  element  geometry.  The  program  can  be  used  to  conduct  plane 
stress,  plane  strain,  or  axisymnetric  analysis  of  an  adhesively 
bonded  structure  subject  to  a  time  varying  thermal  and  mechanical 
loading.  A  nonlinear  Fickean  moisture  diffusion  model  and  a  energy- 
based  delayed  failure  criterion  are  also  provided  in  NOVA. 


1.3.  A  Review  of  the  Literature 

Adhesive  bonding  has  been  in  use  for  many  years.  Most  of  the 
early  bonds  used  animal  and  vegetable  glues,  and  the  structural  use 
of  these  glues  has  been  confined  mostly  to  timber.  The  use  of 
synthetic  resins  in  the  structural  bonding  of  timber  began  in  early 
1930' s.  Synthetic  resins  are  less  susceptible  to  moisture,  fungi  and 
bacteria.  In  recent  years,  synthetic  polymers,  because  of  their 
strength  and  other  bonding  properties  have  been  widely  used  to  bond 
metals  and  composite  materials.  As  noted  in  the  introduction, 
adhesive  joints  have  several  advantages  over  the  conventional 
mechanical  fasteners  (e.g.  bolts,  rivets  and  welding).  These 
include:  lighter  weight,  savings  in  production  cost,  avoid  stress 
concentrations  and  thermal  distortions  due  to  the  hole  drilling  or 
welding,  and  bonding  of  dissimilar  and/or  brittle  materials. 
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An  analysis  of  adhesive  stresses  in  bonded  joints  which  Included 
the  effects  of  load  eccentricity  was  first  performed  by  Go land  and 
Relssner  13]  In  1944  under  the  following  assumptions: 

1.  Adhesive  flexibility  Is  negligible,  and  joint  is  homogeneous 
(i.e..  Ignore  the  presence  of  the  adhesive), 

2.  No  axial  stress  exists,  and  other  stresses  do  not  vary 
through  the  thickness  of  the  adhesive  layer. 

Under  these  simplifying  assumptions,  Goland  and  Relssner  [3] 
developed  one-dimensional  elasticity  solutions  for  two  limiting 
cases:  (1)  the  case  In  which  the  adhesive  layer  Is  homogeneous,  thin 
and  stiff  so  that  its  deformation  can  be  neglected,  the  axial  stress 
Is  zero  and  stresses  do  not  vary  through  the  adhesive  layer;  and  (11) 
the  case  In  which  the  adhesive  layer  Is  soft  and  flexible  and  the 
joint  flexibility  Is  mainly  due  to  the  deformation  of  the  adhesive 
layer  (I.e.,  adherends  are  rigid),  the  axial  stress  is  zero  and 
stresses  do  not  vary  through  the  adhesive  layer.  In  the  first  case, 
the  peel  stress  is  found  to  be  very  high  at  the  edge  of  the  joint, 
while  the  shear  stress  Is  zero.  In  the  second  case,  the  maximum 
values  of  the  peel  and  shear  stresses  occur  at  the  edges  of  the 
joint.  The  Goland-Reissner  analysis  gives  resonable  solutions  for 
shear  and  peel  stresses  in  the  interior  of  the  adhesive  layers;  the 
Joint-edge  loads  are  not  in  equilibrium. 

Erdogan  and  Ratwani  (4|  presented  analytical  solution  based  on  a 


\ 


one-dmensional  model  for  calculating  stresses  in  a  stepped  lap 


Joint.  One  adherend  was  treated  as  Isotropic  and  the  second  as 
orthotropic,  and  linear  elastic  behavior  was  assumed.  The  thickness 
variation  of  the  stresses  in  both  the  adherends  and  in  the  adhesive 
was  neglected. 

Wooley  and  Carver  (51  determined  rtress  distributions  in  a 
simple  lap  joint  using  the  finite  element  method.  The  constant 
strain  quadrilateral  obtained  by  combining  four  constant  strain 
triangular  elements  was  used.  One  end  of  the  adherend  was  assumed  to 
be  hinged  and  other  end  was  allowed  to  move  freely  In  the  direction 
parallel  to  the  original  bond  line.  They  Investigated  the  influence 
of  the  ratio  of  the  Young's  moduli  of  adherend  to  adhesive  materials 
and  geometries  on  the  peel  and  shear  stress  distributions.  The 
results  compared  favorably  with  the  results  of  Goland  and  Reissner. 

Hart-Smith  (6)  Improved  upon  the  approach  of  Goland  and  Reissner 
by  considering  a  third  free-body-dlagram  for  the  adherend  outside  the 
joint  in  addition  to  the  two  free-body-diagrams  from  each  of  the 
upper  and  lower  halves  of  the  joint.  With  three  separate  sections  to 
consider,  three  relations  between  displacements  and  bending  moments 
were  obtained.  Additional  boundary  conditions  involving 
displacements  and  their  first  derivatives,  not  considered  by  Goland 
and  Reissner,  were  imposed  in  order  to  solve  for  the  additional 
unknowns.  In  addition  to  the  improvement  on  the  analysis  of  Goland 
and  Reissner,  Hart-Smith  (61  also  established  a  quantitative 
influence  of  adhesive  plasticity  in  shear.  The  elastic-plastic 
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theory  used  by  Hart-Stn1th  predicts  an  increase  in  joint  strength,  and 
was  shown  to  be  capable  of  explaining  premature  failure  predictions 
found  when  using  linear  elastic  analyses.  The  quantitative  effects 
of  stiffness  imbalance  were  also  accounted  for. 

A  finite-element  stress  analysis  for  adhesive  lap  joints  using 
linear  elasticity  and  elasto-plasticity  theories  was  reported  by  Liu 
(7).  Stress  distributions  in  the  adhesive  layer  for  different  joint 
parameters  (geometry,  material  properties,  and  loading  conditions) 
were  studied  and  compared.  The  existence  of  stress  gradients  through 
the  thickness  of  the  adhesive  layer,  close  to  the  joint  edges,  was 
observed  by  Adams  and  Peppiatt  [81.  They  subsequently  performed  a 
linear  elastic  finite  element  analysis  on  adhesively  bonded  lap 
joints  employing  more  than  one  element  through  the  thickness  of  the 
adhesive  layer,  close  to  the  joint  edges.  Adams  and  Peppiatt  [9] 
also  studied  the  adhesive  yielding  in  double  bevel  and  scarf 
joints.  The  adhesive  was  assumed  to  be  elastic-perfectly  plastic.  A 
nonlinear  finite-element  analysis  of  single  and  double  lap  joints  was 
presented  by  Humphreys  and  Herakovich  110].  The  nonlinear  stress- 
strain  response  was  represented  by  the  Ramberg-Osgood 
approximation.  Mechanical  and  thermal  loadings  were  considered  but 
only  one  element  through  the  thickness  of  the  adhesive  layer  was 
used. 

Allman  [11)  presented  an  elastic  stress  analysis  based  on  the 
strain  energy  density  of  a  particular  joint.  The  effects  of  bending. 
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stretching  and  shearing  of  the  adherends  were  included,  and  the 
shearing  and  tearing  action  In  the  adhesive  was  accounted  for.  All 
conditions  of  stress  equilibrium  In  the  joint  and  stress-free  surface 
conditions  were  satisfied.  It  was  assumed,  however,  that  the  axial 
stress  varies  linearly  through  the  adherend  thicknesses  and  that  the 
shear  stress  Is  constant  through  the  adhesive  thickness.  Allman 
obtained  solutions  for  the  single  lap  joint,  although  the  method  also 
appears  to  be  applicable  to  other  joint  configurations.  He  found 
that  the  average  shear  stress  concentration  is  11*  higher  than  that 
of  Goland  and  Relssner's  first  analysis,  while  the  average  peel 
stress  at  the  joint  edge  is  67*  lower.  Compared  with  the  second 
analysis  of  Goland  and  Relssner,  Allman's  method  yielded  a  shear 
stress  concentration  of  15*  and  31*  less  for  metal  and  composite 
adherends,  respectively,  while  the  average  peel  stress  at  the  joint 
edge  was  27*  higher  and  36*  lower  for  the  same  types  of  adherends, 
respectively. 

Phenomenological  considerations  were  discussed  by  Hart-Smith 
[121  which  greatly  Improve  our  understanding  of  the  sources  of  non- 
uniform  load  transfer,  viz.,  adherend  extensivity,  stiffness 
imbalance  and  thermal  mismatch.  He  also  explained  how  the  lightly 
loaded  central  area  of  the  joint,  away  from  the  joint  edges, 
restricts  cumulative  creep  damage,  and  suggests  that  this  region  is 
vital  for  long  term  durability.  The  amount  of  lightly  loaded  central 
area  is  a  function  of  the  overlap  length. 
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Yuceoglu  and  Updike  [13]  presented  a  numerical  method  for 
solving  peel  and  shear  stresses  in  the  adhesive  of  double  lap,  double 
strap  and  stiffner  plate  joints.  Bending  and  transverse  shear  were 
included  in  the  analytical  model.  Shear  stresses  were  not  required 
to  drop  to  zero  at  the  joint  edges  after  reaching  peak  values  close 
to  the  edges.  Yuceoglu  and  Updike  maintained  that  an  analytical 
model  which  would  allow  the  shear  stresses  to  drop  to  zero  at  the 
joint  edges  would  give  approximately  the  same  or  slightly  lower  peak 
values  of  shear  and  peel  stresses.  Their  method  also  reveals  that 
adherend  bending  has  a  significant  effect  on  both  adhesive  shear  and 
peel  stresses,  especially  the  latter. 

Delale  and  Erdogan  [14,15]  performed  a  plate  analysis  similar  to 
that  of  Goland  and  Reissner  on  the  single  lap  joint  assuming  linear 
elastic  adherends  and  a  linear  viscoelastic  adhesive.  Separate 
stress  distributions  were  calculated  for  membrane  loading,  bending, 
and  transverse  shear  loading.  They  further  extended  their 
viscoelastic  analysis  of  the  single  lap  joint  to  include  time- 
dependent  temperature  variations  approximated  by  a  piecewise  constant 
function. 

6ali  and  Ishai  [16]  perfornwd  a  finite  element  analysis  on  a 
symmetric  doubler  model  with  linear  elastic  adherends  and  the 
adhesive  obeying  a  nonlinear  effective-stress-strain  relationship. 

The  effective-stress-strain  relationship  was  derived  from  stress- 
strain  curves  obtained  by  tensile  and  shear  test  data,  and  based  on 


the  Von  Mises  deviatoric  energy  yield  criterion.  An  iteration 
procedure  was  applied  to  the  linearly  elastic  finite  element  problem 
using  a  specific  secant  modulus  for  each  element  separately.  The 
secant  modulus  was  found  from  the  corresponding  effective  strain  of 
the  previous  solution  and  the  corresponding  effective  stress  was 
found  from  the  experimental  stress-strain  curves.  Gall  and  Ishai 
analyzed  their  symmetric  doubler  model  using  both  plane  stress  and 
plane  strain  and  found  that  the  plane  strain  solutions  converged 
faster  and  yielded  less  conservative  results,  i.e.,  lower  stresses, 
than  the  plane  stress  solutions.  Nonlinear  solutions  were  also  found 
to  be  considerably  lower  than  the  linear  solutions,  the  difference 
being  more  pronounced  in  the  plane  stress  case.  The  problem  was  also 
solved  with  the  adhesive  following  an  elastic-perfectly-plastic 
effective-stress-strain  law.  The  difference  between  these  results 
and  those  of  the  continuous  nonlinear  effective-stress-strain  case 
was  found  to  be  very  smal 1 . 

Nagaraja  and  Alwar  (171  analyzed  a  tubular  lap  joint  with  the 
finite  element  method  assuming  linear  elastic  adherends  and  a 
nonlinear  biaxial  stress-strain  law  in  the  adhesive.  The  constants 
appearing  in  the  nonlinear  law  were  obtained  from  uniaxial  tension 
test  data.  The  stress-strain  relationship,  however,  was  assumed  to 
be  time- independent.  Nagaraja  and  Alwar  demonstrated  that  for  low 
stress  levels,  of  the  order  of  12X  of  the  fracture  stress,  the 
nonlinear  stresses  were  as  much  as  15X  lower  in  shear  and  8X  lower  in 


peel  than  the  linear  stresses.  Nagaraja  and  Alwar  [18]  also 
performed  a  finite  element  analysis  on  a  single  lap  joint,  treating 
the  adherends  as  linear  elastic  materials  but  the  adhesive  as  a 
linear  viscoelastic  material.  The  relaxation  modulus  was  assumed  to 
be  equal  to  the  inverse  of  the  creep  compliance,  the  latter  being 
obtained  experimentally. 

Only  recently,  work  involving  the  time  dependent  fracture 
characteristics  of  adhesively  bonded  joints  has  been  under  way. 
Francis  et  a1.  (19)  discussed  the  effects  of  a  viscoelastic  adhesive 
layer,  geometry,  mixed  mode  fracture  response,  mechanical  load 
history,  environmental  history  and  processing  variations  on  the 
fracture  processes  of  adhesively  bonded  joints.  However,  their 
finite  element  analysis  Includes  only  linear  elastic  fracture 
mechanics. 

Dattaguru,  et  al.  [201  have  performed  cyclic  de-bond  research  on 
the  crack  lap  specimen  and  performed  analyses  with  a  finite  element 
program  GAMMAS,  developed  in-house  at  NASA-Langley.  Their  program 
Includes  geometric  and  material  nonlinearities  but  does  not  include 
viscoelastic  capability.  Also,  fracture  is  modeled  using  linear 
elastic  fracture  mechanics  but  no  failure  law  is  included. 

Botha,  Jones  and  Brinson  (21),  Henriksen  [22],  Becker,  et  al . 
(23),  and  Yadagiri  and  Papi  Reddy  (24)  reported  results  of 
viscoelastic  finite-element  analysis  of  adhesive  joints.  Henriksen 
used  Schapery's  (25)  nonlinear  viscoelastic  model  to  verify  the 


experimental  results  of  Peretz  and  Weitsman  [261  for  an  adhesive 
layer.  The  work  of  Becker  et  al.  (231  is  largely  concerned  with  the 
development  of  a  finite-element  stress  analysis  program,  called 
VISTA,  for  adhesively  bonded  joints.  The  'intrinsic  nonlinear  model' 
based  on  free-volume  concept  of  Knauss  and  Emri  [27)  was  used  in 
VISTA.  The  work  of  Yadagiri  and  Papi  Reddy  (24)  is  limited  to  linear 
viscoelastic  analysis.  Botha  et  al.  [21]  considered  linear  and 
bilinear  viscoelastic  models  in  their  study. 

Pickett  and  Hollaway  [28]  presented  both  classical  and  finite 
element  solutions  for  elastic-plastic  adhesive  stress  distributions 
in  bonded  lap  joints.  Single,  double  and  tubular  lap  configurations 
having  both  similar  and  dissimilar  adherends  were  considered.  The 
results  show  how  the  development  of  adhesive  yielding  will  occur  as 
the  joints  are  loaded  to  a  failure  condition.  The  detrimental  effect 
of  adherend-stiffness-imbalance  on  the  adhesive  stress  distribution 
was  also  shown. 

An  approximate  method  to  analyze  viscoelastic  problems  has  been 
outlined  by  Schapery  [291.  In  this  method,  the  solution  to  a 
viscoelastic  problem  is  approximated  by  a  corresonding  elasticity 
solution  wherein  the  elastic  constants  have  been  replaced  by  time 
dependent  creep  or  relaxation  functions.  The  method  may  be  applied 
to  linear  as  well  as  nonlinear  problems.  Weitsman  [30]  used 
Schapery's  quasi-elastic  approximation  to  investigate  the  effects  of 
nonlinear  viscoelasticity  on  load  transfer  in  a  symmetric  double  lap 


joint.  Introducing  a  stress-dependent  shift  factor,  he  observed  that 
the  enhanced  creep  causes  shear  stress  relief  near  the  edges  of  the 
adhesive  joint. 

Alvazzadeh  and  Verchery  [311  developed  several  special  interface 
finite  elements  based  on  Relssner's  principle  which  take  into  account 
all  the  continuity  conditions  at  the  adhesive-adherend  interface. 
These  elements  were  then  used  to  perform  a  two-dimensional  analysis 
of  an  adhesively  bonded  butt  joint.  It  was  observed  that  the 
interface  stress  distribution  could  be  evaluated  more  accurately 
using  the  interface  finite  elements  compared  to  classical  ones. 

The  analytical  procedures  reviewed  in  the  preceding  pages  are 
primarily  applicable  for  bonded  joints  with  homogeneous  isotropic 
adherends.  These  procedures  have  been  modified  for  composite 
adherends  to  account  for  their  anisotropic  and  heterogeneous 
nature.  A  comprehensive  review  of  publications  relating  to  all 
aspects  of  adhesively  bonded  joints  in  composite  materials  is 
presented  in  [32]. 

Reddy  and  Sinha  [331  extended  the  work  of  Erdogan  and  Ratwani 
[41  to  obtain  analytically,  the  stress  distribution  in  adhesively 
bonded  joints  between  two  orthotropic  materials.  Similarly,  Renton 
and  Vinson  [341  extended  the  work  of  Goland  and  Reissner  [31  to 
obtain  the  linear  elastic  response  of  two  generally  orthotropic 
adherends  adhesively  bonded  together.  Barker  and  Hatt  [35]  used 
linear  elastic  finite-element  analysis  to  evaluate  the  behavior  of  an 
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adhesive  joint  bonding  an  advanced  composite  to  a  metallic 
substrate.  A  special  element  was  used  to  model  the  adhesive  layer  as 
a  separate  elastic  medium  of  finite  thickness  in  order  to  remove  the 
stress  singularity  that  exists  when  dissimilar  materials  are 
joined.  Adams  and  Peppiatt  [36]  performed  two  dimensional  finite 
element  stress  analysis  on  lap,  bevel  and  scarf  joints.  The 
adherends  were  treated  as  homogeneous  anisotropic  materials  with 
linear  elastic  properties,  while  the  adhesive  was  treated  as  an 
elastic-perfectly  plastic  material.  The  effect  of  adhesive  spew 
fillets  on  the  stress  distribution  was  also  taken  into  account  in 
this  analysis.  It  was  observed  that  the  predicted  joint  efficiency 
was  almost  doubled  when  non-linear  adhesive  behavior  was  accounted 
for. 

Renton  and  Vinson  [37]  used  laminated  plate  theory  coupled  with 
a  structural  mechanics  approach  to  obtain  analytical  solutions  for 
stresses  and  deformations  within  a  bonded  single  lap  joint.  The 
closed  form  solutions  were  then  used  to  conduct  a  parametric  study 
which  revealed  that  changes  in  adhesive  moduli,  adherend  longitudinal 
modulus,  and  bond  overlap  length  have  a  significant  effect  on  the 
magnitude  of  the  peak  stresses  within  the  adhesive.  However,  changes 
in  adherend  ply  orientation  appeared  to  have  only  a  modest  influence 
on  adhesive  peak  stresses.  These  predictions  were  verified  by  the 
authors  from  experimental  observations. 
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Schaffer  and  Adams  (381  carried  out  a  nonlinear  viscoelastic 
analysis  of  a  unidirectional  composite  laminate  using  the  finite 
element  method.  The  nonlinear  viscoelastic  constitutive  law  proposed 
by  Schapery  [251  was  used  in  conjunction  with  elastoplastic 
constitutive  relations  to  model  the  composite  response  beyond  the 
elastic  limit. 

Ghoneim  and  Chen  (391  developed  a  viscoelastic-viscoplastic  law 
based  on  the  assumption  that  the  total  strain  rate  tensor  can  be 
decomposed  into  a  viscoelastic  and  a  viscoplastic  component.  A 
linear  viscoelasticity  model  is  used  in  conjunction  with  a  modified 
plasticity  model  in  which  hardening  is  assumed  to  be  a  function  of 
viscoplastic  strains  as  well  as  the  total  strain  rate.  The  resulting 
finite  element  algorithm  is  then  used  to  analyze  the  strain  rate  and 
pressure  effects  on  the  mechanical  behavior  of  a  viscoelastic- 
viscoplastic  material. 

Analysis  of  crack  growth  in  viscoelastic  media  are  mainly 
limited  to  linear  isotropic,  homogeneous  materials.  Schapery  [40] 
proposed  the  use  of  parameters  similar  to  the  J  integral  for  quasi¬ 
static  crack  growth  in  a  class  of  nonlinear  viscoelastic  materials 
subject  to  finite  strains. 

Czarnocki  and  Piekarski  (41)  used  a  nonlinear  elastic  stress- 
strain  law  for  three-dimensional  failure  analysis  of  a  symmetric  lap 
joint.  Taking  into  account  the  variation  of  Poisson's  ratio  with 
strain  within  the  adhesive,  the  authors  concluded  that  the  failure  of 
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the  adhesive  layer  originates  in  the  central  plane  of  a  joint  (at  the 
front  edge).  It  was  also  observed  that  the  joint  width  does  not  have 
any  effect  on  the  stress  peaks  in  the  central  plane  and  that  the 
application  of  a  weaker  but  more  flexible  adhesive  results  in  higher 
load  carrying  capacity  and  lower  stress  concentrations  in  the 
adherends. 

Over  the  years  several  time-dependent  failure  criteria  have  been 
proposed  for  predicting  yield  and  fracture  of  polymeric  materials. 
Nagdhi  and  Murch  (421  and  Crochet  (43)  have  used  a  modified  von  Mises 
criteria  for  viscoelastic  materials  by  assuming  that  the  radius  of 
the  yield  surface  depends  upon  the  strain  history.  An  energy  based 
delayed  failure  criterion  for  polymeric  materials  subjected  to 
isothermal  creep  was  developed  by  Reiner  and  Weissenberg  (44). 
According  to  this  theory,  failure  occurs  when  the  stored  deviatoric 
strain  energy  in  a  viscoelastic  material  reaches  a  certain  maximum 
value  called  the  resilience,  which  is  a  material  constant.  Bruller 
(45)  and  Hiel  et  al.  (46)  applied  the  Reiner-Weisenberg  failure 
criterion  to  various  viscoelastic  materials,  including  composite 
laminates,  and  obtained  good  agreement  with  experimental 
observations. 

It  is  now  well  known  that  moisture  diffusion  can  have  a 
significant  effect  on  the  stress  field  within  an  adhesive  layer  in  a 
bond.  Weitsman  (47)  used  a  variational  method  coupled  with  Fickean 
diffusion  law  to  study  the  interfacial  stresses  in  viscoelastic 
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adhesive- layers  due  to  niolsture  sorption.  From  the  results  of  this 
uncoupled  linear  thermoelastic  analysis,  he  concluded  that  the 
location  of  the  maximum  interfacial  tensile  stress  depends  on  the 
geometry  of  the  joint  as  well  as  the  progress  of  the  diffusion 
process  within  the  joint.  Weitsman  (48)  used  the  correspondence 
principle  to  generate  a  linear  viscoelastic  solution  from  the  linear 
elastic  analysis  of  moisture  sorption  within  an  adhesive  layer.  He 
observes  that  the  viscoelastic  analysis  predicts  detrimental  effects 
due  to  stress  reversals  caused  by  fluctuations  in  relative  humidity, 
that  are  not  predicted  in  an  elastic  analysis.  However,  he 
acknowledged  the  omission  of  the  effect  of  moisture  content  on  the 
viscoelastic  response  of  the  resins  in  his  analysis. 

Tobing,  et  a1.  (491  used  the  finite  element  method  to  study  the 
micro-mechanical  effect  of  moisture  sorption  in  graphite-epoxy 
composites.  Using  a  constitutive  equation  based  on  the  Flory-Huggins 
lattice  model  for  polymer  solvent  interactions,  they  concluded  that 
the  stresses  at  the  graphite-epoxy  interface  have  a  strong  dependence 
on  moisture  content,  fiber  spacing,  and  applied  load. 

Yaniv  and  Ishai  (50)  developed  a  linear  viscoelastic  closed  form 
solution  as  well  as  a  nonlinear  finite  element  solution  algorithm  to 
study  the  hygrothermal  effects  in  a  bonded  fiber-reinforced 
plastic/aluminum  system.  The  linear  solution  was  used  for  short-term 
predictions  at  low  strain  levels,  whereas  the  finite  element  solution 
was  used  for  long  term  predictions  in  which  geometric  and  material 
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nonlinearities  were  taken  into  account.  The  authors  observed  that 
the  presence  of  moisture  tends  to  considerably  reduce  the  stress 
level  in  the  adhesive  layer  and  may  lead  to  significant  variation  in 
the  time-dependent  deformation  of  the  test  specimen  as  compared  to 
the  dry  state. 

In  the  references  cited  above,  various  authors  have  underscored 
the  effect  of  moisture  content  on  the  viscoelastic  response  of  the 
test  specimen.  However,  the  effect  of  the  viscoelastic  stress  field 
on  the  diffusion  coefficient  was  not  considered.  Lefebvre  et  al. 

[51]  extended  the  free  volume  concept  to  define  a  diffusion 
coefficient  that  is  a  function  of  temperature,  dilatational  strain 
and  solvent  concentration.  The  proposed  nonlinear  diffusion  model 
showed  good  predictive  capability  for  different  values  of  temperature 
and  moisture  conentrations.  They  concluded  that  in  order  to  obtain 
an  accurate  solution  for  the  hygrothermal  effects  within  an  adhesive 
bond,  the  nonlinear  diffusion  problem  needs  to  be  solved  in 
conjunction  with  the  nonlinear  viscoelastic  boundary-value  problem 
until  convergence  is  achieved. 

A  review  of  the  literature  reveals  that  previous  finite-element 
analyses  of  adhesive  joints  were  either  based  on  simplified 
theoretical  models  or  the  analyses  themselves  did  not  exploit  the 
full  potential  of  the  finite  element  method.  Also,  several 
investigations  involving  finite  element  analyses  of  the  same  adhesive 
joint  have  reported  apparent  contradicting  conclusions  about  the 
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variations  of  stresses  in  the  joint  (24,521.  Thus,  there  is  a  need 
for  a  closer  examination  of  the  theories,  underlying  assumptions  on 
material  behavior  and  boundary  conditions,  and  the  finite  element 
formulations  used  in  the  analyses  of  adhesive  joints. 


CHAPTER  2 


NONLINEAR  DESCRIPTION  OF  SOLIDS 
2.1  Introduction 

In  the  linear  description  of  the  motion  of  solid  bodies  it  is 
assumed  that  the  displacements  and  their  gradients  are  infinitely 
small  and  that  the  material  is  linearly  elastic.  In  addition,  it  is 
also  assumed  that  the  nature  of  the  boundary  conditions  remains 
unchanged  during  the  entire  deformation  process.  These  assumptions 
imply  that  the  displacement  vector  u  Is  a  linear  function  of  the 
applied  load  vector  F  ,  i.e.,  if  the  applied  load  vector  is  a  scalar 
multiple  aF  then  the  corresponding  displacements  are  au  . 

The  nonlinearity  in  solid  mechanics  arises  from  two  distinct 
sources.  One  due  to  the  kinematics  of  deformation  of  the  body  and 
the  other  from  constitutive  behavior  (e.g.,  stress-strain 
relations).  The  analyses  in  which  the  first  type  of  nonlinearity  Is 
considered  are  called  geometrically  nonlinear  analyses,  and  those 
In  which  the  second  type  are  considered  are  called  materially 
nonlinear  analyses.  The  geometrically  nonlinear  analysis  can  be 
subclassified  according  to  type  of  nonlinearities  considered.  Two 
such  cases  are;  (1)  large  displacements,  large  rotations,  but  small 
strains,  and  (11)  large  displacements,  large  rotations  and  large 
strains.  In  the  first  case  it  is  assumed  that  rotations  of  line 
elements  are  large,  but  their  extensions  and  changes  of  angles 
between  two  line  elements  are  small.  In  the  second  case  the 
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extension  of  a  line  element  and  angle  changes  between  two  line 
elements  are  large,  and  displacements  and  rotations  of  a  line  element 
are  also  large. 

2.2  Incremental  Equations  of  Motion 

In  the  Lagrangian  description  of  motion  all  variables  are 
referred  to  a  reference  configuration,  which  can  be  the  initial 
configuration  or  any  other  convenient  configuration.  The  description 
In  which  all  variables  are  referred  to  the  Initial  configuration  Is 
called  the  total  Lagrangian  description,  and  the  one  in  which 
all  variables  are  referred  to  current  configuration  Is  called  the 
updated  Lagrangian  description. 

The  equations  of  the  Lagrangian  Incremental  description  of 
motion  can  be  derived  from  the  principles  of  virtual  work  (i.e., 
virtual  displacements,  virtual  forces  or  mixed  virtual  displacements 
and  forces).  Since  our  ultimate  objective  is  to  develop  the  finite 
element  model  of  the  equations  governing  a  body,  we  will  not  actually 
derive  the  differential  equations  of  motion  but  utilize  the  virtual 
work  statements  to  develop  the  finite  element  models. 

The  displacement  finite  element  model  is  based  on  the 
principle  of  virtual  displacements.  The  principle  requires  that  the 
sum  of  the  external  virtual  work  done  on  a  body  and  the  internal 
virtual  work  stored  in  the  body  should  be  equal  to  zero  (see  [53]): 
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J  4(26^^)  dV  -  4(^F)  .  0  (2.1) 

^2 


where 


»  J  du^  dV  +  ;  6u^  dS 


=  the  Cartesian  components  of  the  Cauchy  stress  tensor  in 
configuration  C2  at  time  (t  +  At)  occupying  the  volume 


=  the  Cartesian  components  of  the  infinitesimal  strain 
tensor  associated  with  the  displacements  u^  in  going 
from  configuration  Cj^  at  time  t  to  C2  at 
time  (t  +  At)  : 


1  3U.  3U. 

a  ss  ~  J  \ 

2®1j  2  ^3Xj  3X^^ 

Cartesian  components  of  a  point  in  configuration 


(2.2) 


f^  =  Cartesian  components  of  the  body  force  vector 
measured  in  C2 

2 

t^  =  Cartesian  components  of  the  surface  stress  vector 
measured  in  C2. 

Here  6  denotes  the  variational  symbol  (i.e.,  6u^  denotes  the  virtual 
displacement  in  u^)  and  dV  and  dS  denote  the  volume  and  surface 
elements  in  configuration  over  which  the  integrals  are  defined. 
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Equation  (2.1)  Is  not  so  useful  in  Its  present  form  because  the 
Integrals  are  defined  over  the  volume  V2  and  surface  $2  of  the 
configuration  which  is  yet  unknown.  In  the  linear  analysis,  it 
Is  assumed  that  the  configuration  of  the  body  remains  unchanged  and 
therefore  Eq.  (2.1)  applies  to  the  initial  (undeformed) 
configuration.  The  fact  that  the  configuration  of  the  body  changes 
continuously  In  a  nonlinear  analysis  requires  us  to  use  appropriate 
measures  of  stress  and  strain  and  their  Interrelationship  (i.e., 
constitutive  equations)  so  that  Eq.  (2.1)  can  be  used  to  calculate 
the  configuration  C^.  In  order  to  compute  the  current  configuration 
C2  (often,  the  displacements,  strains  and  stresses)  from  the 
knowledge  of  applied  forces  and  displacements,  and  known  previous 
configuration  Cj^,  we  must  make  some  assumptions.  A  description  of 
the  procedure  based  on  the  updated  Lagrangian  approach  Is  given 
below. 

The  coordinates  of  a  general  point  In  Cg  and  and  are 
denoted  by  (X°,  X°,  X°);  (X^,  X^,  X3),  and  (x^,  x^.  X3), 
respectively.  The  displacements  of  a  general  point  in  are  denoted 
(^Uj^,  ^u^t  ^03).  In  C2  they  are  given  by 

^u^  =  ^u^  +  u,  ,  i  =  1,2,3  (2.3) 

where  u^  are  the  components  of  the  displacement  vector  from  to  C2. 

During  the  motion  of  the  body,  its  volume,  surface  area, 
density,  stresses  and  strains  change  continuously.  The  stress 
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measure  that  we  shall  use  Is  the  2nd  Piola-KIrchhoff  stress  tensor. 
The  components  of  the  2nd  P1ola-K1rchhoff  stress  tensor  In  Cj^  will  be 
denoted  by  S^j.  To  see  the  meaning  of  the  2nd  Piola-KIrchhoff  stress 
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tensor,  consider  the  force  dF  on  surface  dS  In  C2.  The  Cauchy  stress 
tensor  t  Is  defined  by 

(n  .  t)  dS  »  dF  (2.4a) 

A 

where  n  Is  the  unit  normal  to  dS  In  C2.  Mote  that  the  Cauchy  stress 
Is  the  force  per  deformed  area  (I.e.  measured  In  C^)  and  referred  to 
C2.  The  2nd  Piola-KIrchhoff  stress  tensor  at  time  t  +  at  referred  to 
Cj  is  defined  by 

(%  •  \s)  dS^  .  df^  (2.4b) 

A 

where  n^  denotes  the  unit  normal  to  the  surface  element  dS^  in  Cj^. 

The  force  dF^  Is  related  to  dF  by 
-o  - 


where 


dF  -  J"^  •  dF 
-0  - 


,  3X  T 

r  ■  (h>  • 


(2.4c) 


From  the  definition  it  is  clear  that  the  2nd  Piola-Kirchhoff  stress 

is  measured  in  C2  but  referred  to  Cj.  It  can  be  shown  that  that  the 
2  2 

components  and  are  related  according  to 


D  3X* 

0  / _ /  j \ 

0  '3x_'  ^mn  ^3x  ' 
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(2.5a) 
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3X.  -  3X. 
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0  Iff  n 


(2.5b) 


where  denotes  the  density  In  Cj  and  p  Is  the  density  In  C2.  The 

2nd  Plola-KIrchhoff  stress  tensor  Is  syffinetrlc  whenever  the  Cauchy 

2  2  2 

stress  tensor  Is  syffmetrlc.  Note  that  =  2^^j- 

Similarly,  the  Green- Lagrange  strain  tensor  E^j  and  the 
Infinitesimal  strain  tensor  e^^  are  related  by 


‘F  3 _ 3! _ H  p 

rij  3X^  3Xj  2®mn‘ 


(2.6) 


It  Is  also  Important  to  note  that  the  2nd  Plola-KIrchhoff  stress 
tensor  Is  energetically  conjugate  to  the  Green-Lagrange  strain  tensor 
and  the  Cauchy  stress  Is  energetically  conjugate  to  the  infinitesimal 
strain  tensor.  In  other  words,  we  have 


I  .(je,  .)dV. 


Substituting  Eg.  (2.7)  Into  Eg.  (2.1),  we  obtain 
0  -  I  «{^E^j)  dV  -  5(2f). 


(2.7) 


(2.8) 


Next  we  use  the  Incremental  decompositions  of  the  stress  and  strains; 

2c  1  *  c 

l^lj  “  ^1J  l^lj 


i^ij  "  l®ij  ^  Tij 


(2.9) 
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where 


l^U 


I'lj 


incremental  components  of  2nd  Plola-Kirchhoff  stress 
tensor 

(incremental)  components  of  the  infinitesimal  strain 
tensor 
,  au.  au. 

A  1  ^  _ il\ 

2  Vax  )• 


j 


aX 


i 


rij 


1  aX,  aX.* 

*  w 


(2.10) 


Recall  that  is  the  i-th  displacement  component  of  a  generic  point 
in  Cj  (in  going  from  to  C2).  Substituting  Eq.  (2.9)  into  Eq. 
(2.8),  we  have 

“  ■  L  '*'u  ‘‘I'lj  ■"  i^ij' 


or 


-I  K.jSi^e,j)6y.6(^F), 

Linearize  the  equations  by  assuming  that 


l^ij  “  l^ijrs  ®rs*  *  l^ij  “  h®ij 


(2.11) 


(2.12) 


where  are  the  components  of  the  linear  elasticity  tensor. 

Constitutive  equations  of  linear  and  nonlinear  viscoelasticity  will  be 
presented  later.  We  obtain  the  approximate  equation  of  equilibrium. 
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(2.13) 


This  linearization  can  be  Interpreted  as  a  representation  of  the 
nonlinear  curve  between  two  consecutive  load  steps  by  linear  line 
segments. 

2.3  Finite  Element  Model 

Here  we  construct  the  finite  element  model  of  Eq.  (2.13)  for  the 
two-dimensional  case  (see  Reddy  (54)).  Let  each  displacement  increment 
be  approximated  as 

n  . 
aui  »  T  &uJ 

j*  X 

Substituting  Eq.  (2.14)  Into  Eq.  (2.13)  we  obtain 


where 


((<*■]  +  (K*^!)  {AU}  =  {f'-}  -  {F*^} 


H  J  (B^l^  (C)  (B*"!  dA  ,  h  =  thickness 
A, 


^11  ^12  ° 
fCl  =  Cj2  ° 
0  0 


(2.15) 
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^2.2 
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[K°l  =  h  J  (8*^1'  [t]  [B°1  dA 


(t1  * 


Til  ^2  °  ° 

ti2  T22  0  0 

0  0  Til  M2 

0  0  Ti2  T22 


[BM 

(4  X  2n) 


n.l  ° 


{F^}  =  h  J  ('fl'’’  (f}dA.  {F°}  =  h  J  [b'-]’’’  {T}dA 


(4.1 

(2  X  2n) 


il*,  0  III-  0  ...  i|^  0 

I  c  n 

0  i|).  0  41-  ...  0  41 

1  Z 


{t}  =  ^22  * 

M2  ^  ^ 


(2.16) 


It  should  be  noted  that  ^t.^  should  be  computed  using  the  Almansi 


strain  tensor. 
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where. 
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Also,  since  Eq.  (2.13)  is  a  linearized  version  of  Eq.  (2.11),  the  error 
introduced  into  the  calculation  of  the  displacements  u^  between 
configurations  can  drift  the  solution  away  from  the  true  solution 
(especially,  if  the  load  steps  are  large).  Therefore,  a  correction 
should  be  made  to  the  displacements  at  each  load  step.  This  can  be  done 
as  follows:  The  solution  {au>  of  Eq.  (2.15)  allows  us  [with  the  aid  of 
Eq.  (2.3)1  to  compute  the  total  displacements  at  time  (t  +  At)  , 

+  AU^ 

which  can  be  used  to  compute  the  strains  and  stresses  (with  appropriate 
constitutive  equation)  at  time  t  +  At.  By  the  principle  of  virtual 
displacements,  the  true  displacements,  strains  and  stresses  at  any  time, 
say  at  time  t  +  At,  are  such  that  the  internal  virtual  work  is  equal  to 
the  external  virtual  work  done.  Since  au^  (hence  the  strains  and 
stresses  computed  from  them)  are  approximations,  there  will  be  imbalance 
between  the  Internal  and  external  virtual  works  performed  on  the  body. 
This  imbalance  can  be  minimized  by  updating  the  internal  virtual  work 
through  an  iteration  (for  a  fixed  system  of  loads  and  time);  the 
Iteration  is  continued  until  the  imbalance  is  reduced  to  a  preassigned 
value  (i.e.,  a  convergence  limit).  For  example,  displacement  increment 


I 


at  the  (r  +  l)st  Iteration  Is  calculated  from  the  equations 


([K^l  +  {Au}j.  ^  j  =  {fS  - 


(2.18) 


wherein  [t]  and  {t}  are  calculated  using  the  displacements. 


(2.19) 


Equations  (2.18)  and  (2.19)  correspond  to  the  Newton-Raphson 
Iteration.  If  the  left  hand  side  (I.e.  (K°l)  Is  not  updated  during  the 
Iteration,  the  Iterative  scheme  Is  known  as  the  modified  Newton-Raphson 
Iteration. 
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CHAPTER  3 

NONLINEAR  VISCOEUSTIC  FORMULATION 
3.1  Introduction 

A  thermodynamicany  consistent  theory  for  a  single  integral 
representation  of  nonlinear  viscoelasticity  was  first  proposed  by 
Schapery  (25).  The  law  can  be  derived  from  fundamental  principles  using 
the  concepts  of  irreversible  thermodynamics.  A  comprehensive  review  of 
the  thermodynamics  basis  of  Schapery 's  theory  has  been  presented  by  Hiel 
et  al.  [55]. 

The  following  two  sections  deal  with  the  review  and  application  of 
Schapery 's  single  integral  constitutive  law  to  problems  with  uniaxial 
and  multiaxial  states  of  stress  respectively.  The  constitutive 
equations  thus  obtained  are  suitable  for  non-linear  viscoelastic  finite 
element  analysis. 


3.2  Uniaxial  Stress  State 

The  uniaxial  nonlinear  viscoelastic  constitutive  equation  of 
Schapery  [25!  can  be  written  for  an  isotropic  material  as, 


aD(iJi^  -  ij»^)  ^  [g^o^jds 


(3.1) 


In  Eq.  (3.1),  represents  uniaxial  kinematic  strain  at  current  time 
t,  is  the  Cauchy  stress  at  time  t,  0^  is  the  instantaneous  elastic 
compliance  and  &D(\j»)  is  a  transient  creep  compliance  function. 
Superscript,  t,  denotes  current  time.  The  factor  g^  defines  stress  and 
temperature  effects  on  the  instantaneous  elastic  compliance  and  is  a 


m 
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measure  of  state  dependent  reduction  (or  increase)  in  stiffness, 

^0  “  ^’^®"sient  (or  creep)  compliance  factor  gj  has  similar 

meaning,  operating  on  the  creep  compliance  component.  The  factor 
g2  accounts  for  the  influence  of  load  rate  on  creep,  and  depends  on 
stress  and  temperature.  The  function  represents  a  reduced  time  scale 
parameter  defined  by. 


(3.2) 


where,  a^y  is  a  time  scale  'shift  factor'.  For  thermorheologically 
simple  materials,  a  »  a(T)  is  a  function  of  temperature  T  only.  This 
function  modifies,  in  general,  viscoelastic  response  as  a  function  of 
temperature  and  stress.  Mathematically,  a^y  shifts  the  creep  data 
parallel  to  the  time  axis  relative  to  a  master  curve  for  creep  strain 
versus  time.  In  this  model,  four  material  parameters  (g^,  g^,  g^  and  a^) 
are  available  to  characterize  nonlinear  behavior  instead  of  only  one 
with  the  time  scale  shifting  procedure  of  Knauss  and  Emri  [27]. 

The  transient  creep  compliance,  aD(i|;),  can  be  expressed  in  the 
following  exponential  form. 


where 

t 

e  = 


-X'r 

aD(4.)  =  I  D  ll  -  e  1 
r 


(3.3) 


and  are  constants.  Substituting  (3.3)  in  (3.1),  gives, 

lfel9Fl9s. 


.  «r-  iru*rv  ir.-  --v  IT^  ir  ^  -  •  -  -  .  •  ,  *  -  -  ,  ■  V  •  *  m  i^.*  -l  •  •  *  «.  •  -  • 


(3.4) 


L  1 1 1  i  I, 


Wff 
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Letting  the  product  g^o^  be  expressed  as  G^  and  simplifying  the 
Integrand  on  the  right  hand  side  of  Eq.  (3.4)  yields, 

■  ajo/  *Al^  h  -4l\  ![  '*  '  #  >'5- 

(3.5) 

The  third  integration  term  on  the  right  hand  side  of  Eq.  (3.5)  is 
now  separated  into  two  parts,  the  first  part  having  limits  from  zero 
to  (t  -  at)  and  the  second  Integral  spanning  only  the  current  load  step, 
i.e.,  from  (t  -  at)  to  t.  Hence, 


t  _-V(»‘  -  •’)  dG®  ..  r*-“  dG* 


r  e 
0 


ds 


ds 


+  J  e 
t-at 


^ds.  (3.6) 


The  first  term  on  the  right  hand  side  of  Eq.  (3.6)  can  be  rewritten  as. 


t-4t  -X^UW)  jjS 

*  ar 


ds 


t-at  -X  -X  (i|3^"^^-id^) 


dG- 


ds 


-x^(<»^-v^-'^^)  t-at  -x^(^^-^^-,d5)  s 
e  J  e 

0 


ot-at 

=  e  q^ 


ds 


(3.7) 


where. 


I 

I 

hi 


k 

c- 


I 


I 

yj. 


i 

I 


rrirriT* 
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.t-At 


.,t  _  ,t  ,t-At 

All*  =  t|)  -  4i 


^t-At 


ds 


The  second  Integral  on  the  right  hand  side  of  Eq.  (3.6)  is  now 
integrated  by  parts.  To  carry  out  the  integration,  it  is  assumed 
that  varies  linearly  over  the  current  time  step  At.  Hence, 


t-At 


ds 

5 


dG'  e 


ds  X. 


t  ^2-S  'V**  -♦  > 

-  ! 

t-At  t-At  dS^  V 


,  ,  t  ,t-At* 

^1.  dG^-^^  ’ 

dt  x^  dt  x^ 

r  r 


-X  Alp 


,1  -  e 
dt 


_  QVa  f  1.  -  c _ 1 

^  i  \  1 • 


(3. ID) 


In  arriving  at  the  second  step,  the  fact  that  is  assumed  to  be 
linear,  hence  its  second  derivative  is  zero,  is  used.  Since  G^  has  been 
assumed  to  be  a  linear  function  of  time  over  the  current  load 
step  At,  we  can  write, 

dG^  G^-G^-^"^ 


dt 


,t  ,t-At 
P  -  Ip 


I 


& 

(3.8) 

i 

■  xJ 

s 

ds. 

(3.9) 

'A 

iiS 

I 

A 

s 


r 

I 


I 


» 


:  vl: 


(3.11) 


Substituting  (3.11)  Into  (3.10),  gives 


t  -X  (4»^-V*)  .pS  V  Af  1 

;  e  ^  ds  =  IG^  -  Gt-At^Lz  e 

t-At  X  AU.^ 

r 

t  -X  (iti^-4»^)  .pS  .4.  ^ 

J  e  ^  ds  =  [G  -  G  ]s 

t-At 


-XpAii) 


(3.12) 


where. 


.t  1  -  e-V^^ 


Substituting  (3.9^  and  (3.12)  back  Into  (3.5)  and 
writing  G^  =  g^o^ 

e  =  g^O^a  +  g^  I  D^g^o 

r  Au.^ 

pt  r  n  r/  r  t-At  ^  .  t  t  „t-At  t-At,„t, 

-  91  I  Dple  dp  +  1920  -  92  o  iBpl- 


(3.13) 


(3.14) 


Collecting  those  terms  in  Eq.  (3.14)  that  are  multiplied  by  current 
stress  yields, 

■  l55'’o  *  9i32  I  °r  -  3192  I  I>r8ll»* 

♦  91(1  (3.15) 

r 

Defining  Instantaneous  compliance  0^  as  the  compliance  term 
multiplying  the  instantaneous  stress  o^,  and  the  remaining  terms  in  Eq. 


(3.15)  as  hereditary  strains  E^,  yields. 


=  Ofa^  ■*•  E^  =  F(o) 


(3.16) 


where. 


°I  =  9^0  ^  9^2  I  °r  -  9^2  I  °r®J  (3.17) 


Hence,  Eq.  (3.16)  expresses  Schapery's  single  integral  constitutive 
law  in  terras  of  a  stress  operator  that  includes  instantaneous  compliance 
and  hereditary  strains. 

It  is  to  be  noted  that  the  terra  in  Eq.  (3.18)  is  the  r^^ 

component  of  the  hereditary  integral  series  at  the  end  of  the  previous 
load  step  (i.e.  at  time  equals  t  -  At).  The  expression  for  the 
hereditary  integral  at  the  end  of  the  current  load  step  (i.e.  at  time  t) 
can  be  derived  in  the  form  of  a  recurrence  formula  as  shown  below. 

By  definition  [see.  Eq.  (3.9)1, 

„t .  r'  dG*  „ 

V  ■  -fg  '  dT 

t-At  -X^(*i^-4»^)  .-S  t  -S 

=  ;  f-ds-J  e  ^^ds.  (3.19) 

0  t-At  9S 

Using  the  results  from  Eqs.  (3.9)  and  (3.12),  the  above  equation  reduces 


where  1s  defined  by  Eq.  (3.13). 


3.3  Multiaxial  Stress  State 

For  a  thermorheologlcal ly  simple  anisotropic  viscoelastic  material 
under  a  multiaxial  state  of  stress,  the  constitutive  law  proposed  by 


Schapery  (251  is. 


3o  . 

—L.  +  -jna  Ae 

3  cf  ^  j  3  (3  i  j  ftin 


(3.21) 


-  »’>  It  <»i  A)'''  (5-«) 

where,  e.jj  and  are  the  strain  and  stress  tensors  respectively,  Gr  is 

*  i  i 

the  Gibbs  free  energy,  o.  .  and  AS„r  are  second  and  fourth  order  material 

^  1  j  mn 

property  tensors  respectively  and  aQ  is  a  material  kernel  function 
defined  in  [25].  The  quantities  G[^,  aQ  and  are,  in  general, 
functions  of  ten  variables,  a.^  and  temperature  T.  Note  that  all 
repeated  indices  in  Eqs.  (3.21)  and  (3.22)  are  to  be  summed  out  over 
their  range  (1,2,3). 

Due  to  the  complex  nature  of  Eq.  (3.21)  it  is  not  possible  to 
determine  the  material  properties  in  this  equation  from  the  uniaxial 
tests  outlined  in  (25).  However,  it  is  possible  to  construct  a  set  of 
small  strain,  three  dimensional  constitutive  equations  from  (3.21), 
which  is  consistent  with  the  thermodynamic  theory  in  [251,  and  which  yet 
enables  all  properties  to  be  evaluated  from  uniaxial  tests.  The 


assumptions  which  need  to  be  made  for  this  purpose  are  as  follows: 


(a)  The  Gibb's  free  energy  G|^  is  a  quadratic  function  of  stress. 


(b)  .,j  =  o,j 

When  the  free  energy  Gp  is  a  quadratic  function  of  stress, 

sGq 

— —  =  sT'Jcoio 

30ij  *  mn 


(3.23) 


(3.24) 


where  S'^j(O)  are  the  instantaneous  components  of  the  linear  viscoelastic 

creep  compliance  tensor.  Equation  (3.24)  implies  that  the  initial 

response  of  the  material  is  linearly  elastic  under  suddenly  applied 

stresses,  which  is  often  the  case  for  metals  and  plastics. 

The  second  assumption,  on  the  other  hand,  leads  to  the 

linearization  of  the  coefficient  of  the  transient  term  in  Eq.  (3.21). 

Mathematically,  this  is  given  by, 

30-.-  ^  ^  and  j  =  n 

^  =  {  (3.25) 

^^ij  0,  if  i  *  m  and  j  *  n 

Equation  (3.25)  implies  that  the  jump  in  strain  due  to  load  application 
equals  the  jump  when  the  load  is  removed.  This  behavior  is  exhibited  by 
some  types  of  plastics  [56]. 

Substituting  Eqs.  (3.24)  and  (3.25)  in  (3.21)  and  (3.22), 


(3.26) 

“5j  ■  L  (3-37) 

where,  g^  =  Superscript,  t,  denotes  current  time.  Equation 

(3.26)  is  a  set  of  three  dimensional  constitutive  equations  for 


► 


anisostropic  viscoelastic  materials  which  includes  the  nonlinear 
functions  and  a^j  appearing  in  the  uniaxial  relations  (3.1)  and 
(3.2).  Note  that  the  functions  g^  and  a^y  are  expressed  as  a  function 
of  the  octrahedral  shear  stress,  which  is  a  stress  invariant. 

For  a  homogeneous  isotropic  nonlinear  viscoelastic  material,  Eq. 
(3.25)  reduces  to  the  form  presented  by  Schapery  ([25]  and  [56]), 

®(J  ■  ♦  (D  -  (3-28) 

where, 

-  *')  It  (3*29) 

{0  -  J}{g2i}  ’  (0(0)  -  0(0)  (-^D(^^  -  i,^) 

.  40(4,^  .  J»")]  (92^0^ 

in  which, 

f  ^  1 

/  -  4."  =  ;  (a^j)  ds  (3.31) 

T 

Expanding  Eq.  (3.28)  term  by  term  for  the  strains, 

*11  *  {0}{g2®ii)  +  {0  -  JHg2‘^ll  92‘^22  *  ^2*^33^ 

Arranging  terms, 

eJi  »  {Olfg^oJll  +  {D  -  J}{g20^2}  +  (D  -  J}{g2a^3} 


(3.32) 


Similarly, 


632  “  (D  -  J}{g2°n}  {0}{g242^  -  ‘^^92^33}  (3.33) 

^^2  =  2{J}{g^oj2}  (3.34) 

633  =  (D  -  (0  -  J}{g^a^2J  ^  {D}(g2°33}  (3-35) 


The  transient  components  of  the  creep  and  shear  compliances  can  be 
written  in  the  form  of  Prony  series  as. 


=10(1-  e'^^^ 

r 

(3.36) 

AJ(^)  =  1  J^Il  -  e  M 

(3.37) 

where  and  are  the  reciprocal  of  the  retardation  times  in  creep  and 
shear  respectively.  Also,  let. 


Considering 
Eqs.  (3.36) 


0(0)  =  0^ 

J(0)  = 

a  term  of  the  form  (Dlfg^o^j.} 
to  (3.39)  gives. 


in  Eq. 


(3.38) 

(3.39) 

(3.32)  and  substituting 


iDKg^j) 


(3.40) 


Recognizing  that  Eq.  (3.40)  is  similar  in  appearance  to  Eq.  (3.6)  and 
making  use  of  the  results  derived  in  Section  3.2, 


(0Kg^«5j)  .  -o-j 


(3.41) 


where  0^ 


is  the  instantaneous  creep  compliance  function  at  time  t,  given 


Of  =  Oq  92  I  -  sp. 


(3.42) 


Q^j  are  the  hereditary  strain  components  due  to  tensile  creep  at  time 


Q,j  -  I  192  Vl 


t-at„t  t-at  t-i 


^  -  e 


(3.43) 


(3.44) 


and,  are  components  of  the  hereditary  Integral  given  by  the 


recurrence  formula. 


=.  nt-at  .  ,_t^t  t-at  t-at,^t  ,, 

^r,ij  ®  V,ij  -  92  °ij  iBp  (3.45) 


Similarly,  a  term  of  the  form  {JKg^oJj.}  in  Eq.  (3.32)  can  be 


expressed  as. 


(3.46) 


where,  Jj  is  the  instantaneous  shear  compliance  function  at  time  t,  and 


is  given  by. 


jf  =  Jo  *  4  J  ■'r*'  -  (J- 

P^j  are  the  hereditary  strain  components  due  to  shear  at  time  t. 


(3.47) 


pt  _  Y  r  t-at  t  t-at 
Pij  I  192  ^^ij 


(3.48) 
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pt  =  I-  e 


(3.49) 


and,  are  components  of  the  hereditary  integral  given  by  the 

recurrence  formula; 


nt  _  t-At  .  t  t  „t-&t  t-at,^t  ,, 

Pr.ij  “  ®  Pr.ij  ^Vij  -  ^2  °ij  IV 


Substituting  Eq.  (3.41)  and  (3.46)  in  Eqs.  (3.32)  to  (3.35),  and 
dropping  superscripts,  gives. 


where. 


®11  " 

“i’ll 

+  (Dj  - 

22  ■"  - 

J,)o33 

*"11 

(3.51) 

®22  " 

(0,  - 

+  DjO 

22  ■ 

'’l'°33 

*  "22 

(3.52) 

n2 

=  2JjC 

'12  ^12 

(3.53) 

®33  " 

(Dl  - 

Jl)0ii 

+  (Oj 

-  Ji)o22  • 

*■  °I'’33 

.H33 

(3.54) 

"u  = 

=  Qll  + 

Q_-  + 
^22 

^33  ■  '’22 

-'’33 

(3.55) 

“22  ' 

'  ^11  ^ 

Q22  + 

^33  ■  '’ll 

-''33 

(3.56) 

”12  " 

2^12 

(3.57) 

"33  ° 

Qll  > 

Q22  > 

^33  -  '’ll 

'  ’’22' 

(3.58) 

Written  in  form  of  a  matrix  equation,  this  becomes. 


WVWA-' 
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{e}  =  [H\{o}  +  {H}. 


(3,59) 


Note  that  the  left  hand  side  of  Eq.  (3.59)  is  a  vector  containing  the 
algebraic  difference  of  kinematic  strains  {e}  and  dllatational  strains 

(er  =  {(«n  -  e).(ej2  -  0),Yi2.(t33  -  9)}  (3.60) 

While,  {a}  contains  four  components  of  Cauchy  stress. 


=  ^‘^ir°22’°12*°33^ 

and  {H}  is  a  vector  of  hereditary  strains  given  by. 

The  matrix  [N]  is  a  4  x  4  coefficient  matrix  given  by. 


(3.61) 


(3.62) 


[N]  = 


0 

(Dl-Jj) 


0 


0 

0 

2Ji 

0 


(D j-J j) 

(Dl-Ji) 

0 

“i 


(3.63) 


Pre-multiplying  Eq.  (3.59)  by  [N]"^,  an  explicit  expression  for  stresses 
in  terms  of  strains  is  obtained, 

{a}  =  [Ml({e}  -  {H})  (3.64) 


I 


."O 


•^1 

Si 


.  -^1 


where. 


[Ml  =  IN)- 


(3.65) 


Equations  (3.64)  and  (3.65)  provide  a  general  viscoelastic 
constitutive  relation  that  can  be  applied  to  either  plane  stress,  plane 
strain  or  axisymmetric  problems.  For  plane  strain,  the  out-of-plane 
strain  component  e33  is  identically  zero.  The  corresponding  stress 
component,  023,  may  be  obtained  from  Eq.  (3.64)  by  setting  e33  =  0. 

Since  for  the  plane  stress  case,  O33  is  identically  zero,  the 
corresponding  strain  component  e33  can  be  evaluated  from  Eqs.  (3.59)  and 


(3.63)  as. 


e-,,  =  (Dt  -  ‘^99) 


(3.66) 


C33  -  -  P'  11  22'  33  vj.wvi/ 

Note  that  the  use  of  creep  and  shear  compliances  as  material  property 
input  allows  the  Poisson's  ratio  to  be  time-dependent.  Hence,  the 
present  formulation  is  applicable  to  any  thermorheological ly  simple 
isotropic  viscoelastic  material  over  any  length  of  time. 

For  the  special  case  where  the  Poisson's  ratio  is  a  constant  with 


time,  then. 


J(i|))  =  (1  +  v)D(il;) 


(3.67) 


The  matrix  [N]  in  Eq.  (3.59)  takes  the  form. 


0  2(l+v) 


(3.68) 
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If  the  viscoelastic  properties  of  a  material  are  defined  by  its 
bulk  and  shear  compliances  instead  of  the  creep  and  shear  compliances, 
then  the  creep  compliance  D(i(»)  in  Eq.  (3.28)  is  replaced  by  the  bulk  and 
shear  properties.  Using  the  viscoelastic  relationship  betwen  creep, 
bulk  and  shear  compliances  given  by. 


D(it»)  =  ^  M(4»)  +  I  J(i-) 


(3.73) 


and  substituting  in  Eq.  (3.28),  the  matrix  relation  given  by  (3.59)  is 
obtained.  However,  for  this  case  the  matrix  [N]  has  the  form. 


[N1  = 


(^Mi.f  Ji) 


(g 


3 


(g  Ml  -  3  Jj) 


(g  Ml  +  ^  Jl) 


(g  Mi  -  3  Jj) 


0 

'I  ■  3 


ih  Mt  -T  Jj) 


0 

0 

2J 

0 


I 


(|  Mj  -^Ji) 
(|  Mj  -ijj) 


(|  Mj  .f  Jj) 


(3.74) 


reditary 

strains  are. 

■ 

11  "  '^11 

-  »(Qjj  .  Q33) 

(3.69) 

■ 

22  "  ^22 

-  m(Qii  *  (J33) 

(3.7Ci 

Hi2  = 

2(1  +  '>)Q32 

(3.71) 

33  "  ^33 

~  '*^^11  ^  ^22^ 

(3.72) 

» 

and  the  corresponding  hereditary  strains  are. 


^11  ^  9  ^22  ■*■  ^33^  3  ^11  "  3  ^'^22  ‘*33^ 

^22  "  i  ^22  ■*■  ^33^  I  ^22  ■  i  ^’*11  *  ’*33^ 


«12  =  ^^2 


(3.77) 


^33  “  9  ^^11  ^22  ^33^  3  ^33  '  3  ^**11  ■*■  ^22^ 


3.4  Finite  Element  Model 

This  section  describes  the  finite  element  implementation  of  the 
nonlinear  viscoelastic  constitutive  law  presented  in  Sections  3.2  and 
3.3.  Since  viscoelastic  materials  often  undergo  large  displacements 
especially  when  subjected  to  creep  type  of  loading,  the  geometrically 
nonlinear  analysis  described  in  Chapter  2  has  been  incorporated  into  the 
viscoelastic  formulation. 

Invoking  the  principle  of  virtual  work  and  following  the  procedure 
outlined  in  Section  2.2  gives, 

''l  ''l 

=  -  f  +  J  ^f^su^dV  +  J  ^.6u.dS  (3.79) 

''l  ''l  h 

where,  are  the  components  of  the  viscoelasticity  constitutive 

tensor.  The  rest  of  the  quantities  and  the  superscripts  in  Eq.  (3.79) 

have  the  same  definitions  as  in  Section  2.2.  Let  each  displacement 

increment  at  any  time  t,  be  approximated  as. 


»y’.:vvyv  vvfv^'.  'Hi.’rw  wvirmr 


‘i '  ji  ''"i*j("i-’‘2)- 


(3.80) 


Substituting  Eq.  (3.80)  into  Eq.  (3.79)  gives, 

((k'-I  +  (K°)){au}  =  {F^}  -  {F‘'}  (3.81) 

where, 

(k'-I  .  hj  Is'-l’^lMIIB'-ldA,  h  .  thickness  (3.82) 

\ 

and  [Ml  is  the  4x4  viscoelastic  constitutive  matrix  defined  in  Eq. 
(3.65).  Note  that  for  the  nonlinear  viscoelastic  case,  the  'linear' 
stiffness  matrix  (K*-)  is  not  really  linear,  but  has  nonlinearities 
imbedded  in  it  due  to  the  presence  of  the  material  kernel  functions  {g^, 
gi,  g2)  in  the  matrix  [Ml.  The  nonlinear  stiffness  matrix  [K°l  is  the 
direct  result  of  the  geometrically  nonlinear  formulation  and  is  given 


[K'^l  =  h/  [B°r[Tl[B°ldA 


(3.83) 


The  definitions  of  matrices  (B*-l,  [B°l,  {F*"},  {F°}  and  [t1  are  the  same 
as  in  Eq.  (2.16).  The  Cauchy  stress  components,  are  computed  by  using 


the  viscoelastic  relation. 


{o}  =  [Ml({e}  -  {H}) 


(3.84) 


which  has  been  derived  in  Section  3.3.  For  a  geometrically  nonlinear 
analysis,  the  vector  {e}  contains  components  of  the  Almansi  strain 


tensor  given  by. 


1  3U.  3U_  3U„  3U 

p  =1  ( _ !i  + _ !!!  _  _ H _ Q\  _  I  9 

km  2  '3X„  3x.  3x.  3X„'  km 

m  k  km 


(3.85) 


It  is  evident  that  Eq.  (3.81)  contains  two  possible  sources  of 
nonlinearities:  material  nonlinearity  due  to  Schapery's  law  and, 
geometric  nonlinearity  arising  from  the  large  displacement 
formulation.  In  order  to  obtain  a  solution  to  this  nonlinear  equation 
at  any  time  step,  the  Newton-Raphson  iterative  technique  is  used.  The 
incremental  displacement  {au}^  obtained  at  the  end  of  the  r^^  iteration 
is  used  to  update  the  total  displacement  for  the  n^^  time  step, 

"V  =  Vl  ^"‘'r 

The  iterative  procedure  continues  until  a  convergence  criterion  is 
satisfied.  After  that,  the  solution  proceeds  to  the  next  time  step. 
Note  that  for  the  first  time-step,  ~ 

The  complete  solution  procedure  for  each  individual  time  step  is 
presented  in  a  logical  step-wise  fashion  and  can  be  used  directly  for 
programming  purposes: 

1.  At  the  beginning  of  each  time  step,  the  stress  vector 

{a}  from  the  previous  time  step  is  accessed.  Note  that  for 
the  initial  or  starting  time  step,  the  stress  vector 
a(t  -  At)  denotes  the  initial  stress  state  at  t  =  0,  given 
by  {a°}.  Since  it  is  customary  to  assume  a  stress  free  state 
to  exist  at  the  start  of  the  solution,  {o°}  is  usually  set  to 
be  zero. 

2.  Temperature  T  at  time  t  is  computed  from  T  =  f(t)  which  is 
supplied  by  the  user  for  problems  involving  thermal  loads. 


3.  The  parameters  g^,  g^,  ^2*  are  known  functions 

of  temperature  and  stress,  are  evaluated  next,  using  the  stress 
vector  obtained  from  previous  load  step. 

4.  Assuming  a^j  to  be  a  linear  function  of  time  over  the  time 
step  At,  the  average  value  of  shift  factor  is  given  as 

a^y  =  change  in  reduced  time  a\i<^ 

is  computed  as  A\Ii^  =  At/a^j  .  In  order  for  this  assumption 

°  avg 

to  be  valid  At  should  be  made  sufficiently  small. 

5.  Hereditary  integral  {q^}  is  computed  using  the  recursive 
formula  given  by  Eq.  (3.45). 

^%xt^  "  ^^^app^  where  x  is  the  load  factor  that  corresponds 
to  the  time  step  under  consideration. 

7.  The  residual  vector  {R}  is  computed  for  each  element  as. 


e 

8.  The  tangent  stiffness  matrix  [K,]®  =  f  [Bl^[M] [B]dv. 

V 

e 

9.  Incremental  displacement  {au}  =  [Ky]"^{R}. 

10.  Total  displacement  {u}^  =  {u}.  ^  +  {au}^.  where  the  subscript  i 
denotes  the  number  of  iterations. 

11.  The  strains  and  stresses  are  computed  using  the  known 
displacement. 

!|{au.}|| 

I|{U^}|| 


12.  Steps  3  through  12  are  repeated  till 


<  tolerance. 


CHAPTER  4 


MOISTURE  DIFFUSION  AND  DELAYED  FAILURE 
4.1  Governing  Equations  for  Diffusion 

The  nonlinear  Fickean  two-dimensional,  diffusion  model  presented  in 
the  present  study  is  the  one  investigated  by  Lefebvre,  et  al.  [51].  The 
diffusion  model  can  also  be  applied  for  penetrants  other  than  moisture. 

Fick's  law  for  the  two-dimensional  diffusion  of  a  penetrant  within 
an  isotropic  material  is  given  by. 


L.  (0  1£)  +  1_  (D  i£)  =  ig 

3x  '  ax^  ay  '  ay^  at 


(4.1) 


where,  c  is  the  penetrant  concentration,  which  is  a  function  of  position 
and  time,  and  D  is  the  diffusion  coefficient. 

In  order  to  model  the  transport  phenomena  in  polymeric  materials, 
Lefebvre  et  al .  [51]  derived  a  nonlinear  diffusion  coefficient  based  on 
the  concept  of  free  volume. 

According  to  this  theory,  the  diffusion  coefficient  for  a  polymeric 
material  above  its  glass  transition  temperature  is  given  by. 


D  =  ^  Te 


-B{l/f  -  1/f  } 


(4.2) 


where,  D  is  the  diffusion  coefficient,  T  is  the  temperature,  f  is  the 
free  volume  fraction,  and  B  is  a  material  constant.  The  subscript  'o' 
denotes  values  at  the  reference  state.  It  is  then  postulated  that  the 
change  in  fractional  free  volume  is  due  to  changes  in  temperature. 
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penetrant  concentration,  and  the  transient  component  of  the  mechanically 
induced  dilatational  strain.  It  is  further  assumed  that  these  changes 
are  additive,  which  is  similar  to  the  assumptions  made  by  Knauss  and 


Emri  [27],  giving. 


f  =  f^  +  3aaT  +  Aej^i^  +  3yC 


(4.3) 


In  Eg.  (4.3),  o  is  the  linear  coefficient  of  thermal  expansion,  y  is  the 
linear  coefficient  of  expansion  due  to  moisture,  N  is  an  exponent  for 
the  saturated  state,  and  ae|^|^  is  the  transient  component  of  the 
mechanically  induced  dilatational  strain.  The  dilatational  strain  due 


to  the  ambient  stress  state  can  be  written  as, 

'kk  ’  ^kkfO)  *  “kk 

where  e|(|((0)  and  Aej^j^  are  the  instantaneous  and  transient  components  of 
the  mechanically  induced  dilatational  strain  en,.  Hence, 


e,,(0)  -  i  M(0)a,, 


=  1  J’g  It  (92‘^kk)'^" 


(4.5) 


(4.6) 


where  M(4i)  is  the  bulk  compliance  of  the  material.  Combining  Eqs.  (4.4) 


and(4.5) , 


^®kk  "  ®kk  ■  3 


(4.7) 


Substituting  (4.7)  in  (4.3), 


f  =  f^  +  3aAT  +  (e^^  -  \  M(0)o^^}  +  3yC^ 


(4.8) 


^  'fc  %i  .  *• 


Substituting  (4.8)  in  (4.2)  gives, 
0. 


rB 


3(aAT  -K  yc^)  +  (e^^  -  i 


D  =  T  exp  (f-  [ - .  ^  - 1}  (4.9) 

0  0  +  3(aAT  +  yC  )  +  (e^^  -  i  M^c^^) 

From  the  viscoelastic  formulation  presented  in  Chapter  3,  it  is  evident 

that  the  dilatational  strain  e|^|^  is  dependent  on  the  stress  history, 

temperature  and  penetrant  concentration,  that  is, 

®kk  "  ^kk^^kk’^’*^^  (4.10) 

Hence,  the  two  sources  of  nonlinearity  in  Eq.  (4.1)  are  moisture 

concentration  c,  and  dilatational  strain  e|^|j.  Consequently,  in  order  to 
accurately  model  the  penetrant  transport  phenomena,  the  diffusion 
boundary-value  problem  needs  to  be  solved  in  conjunction  with  the 
nonlinear  viscoelasticity  boundary-value  problem  by  using  an  iterative 
procedure.  The  same  solution  procedure  can  also  be  applied  for 
diffusion  in  polymeric  materials  where  the  plasticizing  effect  of  the 
penetrant  may  cause  the  viscoelastic  time-scale  shift  factor  to  be 
concentration  dependent,  that  is, 

a^^^  =  a(a,T,c)  (4.11) 

One  example  of  such  a  shift  factor  definition  can  be  found  in  the  work 
of  Knauss  and  Emri  [271,  where  the  authors  have  used  the  unifying 
concept  of  free  volume  to  define  a  shift  factor  given  by, 

aaT  +  yc  +  6e, 


a(T,c,ekk)  =  exp{-  j-  (f^  +  +  yC  + 


where  5  is  the  coefficient  of  the  dilatation  term.  Note  that  the 
coefficients  a  and  y  in  Eqs.  (4.9)  and  (4.12)  are,  in  general,  functions 


of  T,  c  and  e|^|j  but  have  been  assumed  to  be  constant  for  the  sake  of 
simplicity.  This  assumption  Is  valid  for  temperatures  above  the  glass 
transition  temperature. 


4.2  Finite  Element  Formulation 

Pick's  law  for  two  dimensional  diffusion  in  a  homogeneous  isotropic 


material  is  given  by. 


i_  (0  .  i_  (^)  .  inn 

3x  '  3x^  ay  'ay  at 


(4.13) 


subject  to  the  boundary  conditions. 


0  —  n  +  0-^0  +  q  =  0onr.  ,  t>0 

3x  X  3y  y  1 


(4.14) 


c  =  c  on  r2  .  t  >  0 


with  the  initial  condition. 


c  =  c  in  a  ,  t  =  0 
0 


(4.15) 


(4.16) 


where,  a  is  the  two  dimensional  region  in  which  diffusion  occurs, 
and  r  is  the  boundary  to  this  region. 

The  finite  element  formulation  for  Eq.  (4.13)  incorporating  the 
initial  and  boundary  conditions  (Eqs.  (4.14)  to  (4.16))  was  carried  out 
following  the  variational  procedure  used  by  Reddy  [54].  The  time 
dependent  moisture  concentration  is  approximated  as. 


•  .  •r\  * 


(4.17) 


c(x,y,t)  =  I  ■4>.(x,y)c.(t) 
j=l  ^  ^ 


The  resulting  finite  element  equations  cast  in  a  matrix  form  are  given 

by, 

[M<«)|{c)  .  |K<®>|(cl  =  (4.18) 


where. 


.(e)  = 


i 


=  J  •4».ji)4dxdy 

(4.19) 

Q 

34*4  94*.>  3i|>i  34*4 

- + - -^)dxdy 

3X  3X  3y  3y  ^ 

(4.20) 

=  -  J  4»^qds 

(4.21) 

The  superscript  (e)  is  used  to  denote  that  the  equations  are  valid  over 
each  element.  The  range  of  the  subscripts  i  and  j  is  equal  to  the 
number  of  nodes  per  element. 

The  time  derivative  (cj  is  approximated  using  a  0  family  of 
approximations  given  by, 

("Ui  - 1"), 


at 


n+1 


for  0  <  9  <  1  (4.22) 


where,  n  is  the  time  step.  Using  the  approximations  (4.22)  in  (4.18) 
for  time  tp  and  tp^j^  gives, 

-  I8<®'l(c)„  -  (P<=>!„  .  0  (4.23) 


where. 


[b(®)]  =  [M(®h  -  (1  -  0)At  ,[K(®h 


(4.24) 


(4.25) 


{p(®)}  =  +  (1  -  9){f(®^}^)  (4.26) 

Recognizing  that  a  source  of  nonlinearity  in  the  form  of  the  diffusion 
coefficient  0  is  imbedded  in  the  matrix  the  Mewton-Raphson 

technique  is  employed  to  solve  for  the  concentration  at  each  time 

step.  Note  that  for  n  =  1,  the  value  of  (c)  in  Eq.  (4.23)  is  known  from 
initial  conditions. 

4.3  Delayed  Failure:  Uniaxial  Formulation 

When  a  viscoelastic  material  undergoes  deformation,  only  a  part  of 
the  total  deformation  energy  is  stored,  while  the  rest  of  the  energy  is 
dissipated.  This  behavior  is  unlike  elastic  material  where  all  the 
energy  of  deformation  is  stored  as  strain  energy.  Reiner  and  Weisenberg 
[44]  postulated  that  it  is  this  time-dependent  energy 
storage  capacity  that  is  responsible  for  the  transition  from 
viscoelastic  ’•esponse  to  yield  in  ductile  materials  or  fracture  in 
brittle  ones.  According  to  this  theory,  failure  occurs  when  the  stored 
deviatoric  strain  energy  per  unit  volume  in  a  body  reaches  a  certain 
maximum  va'.e  called  the  resilience,  which  is  a  material  property.  Note 
that  when  -.-ere  is  no  dissipation,  that  is,  when  the  material  is 
elastic,  tnen  Reiner-Weisenberg  criterion  becomes  identical  to  the  von 
Mises  criter ■  ::n. 
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Consider  the  single  Kelvin  element  shown  in  Fig.  1,  subject  to  the 
uniaxial  tensile  load  a(t).  The  total  strain  response  e(t)  due  to  the 
applied  stress  can  be  divided  into  two  components:  the  instantaneous 
respone  e^,  and  the  transient  response  ej^(t).  Hence, 

e(t)  =  +  e^(t)  (4.27) 

For  uniaxial  creep,  the  applied  stress  o(t)  is  given  as, 

a(t)  =  a^H(t)  (4.28) 

where  H(t)  is  the  unit  step  function. 

Substituting  Eq.  (4.28)  into  Schapery's  nonlinear  uniaxial  single 
integral  law  given  by  Eq.(3.1),  and  expressing  the  transient  creep 
compliance  D(.('ii)  as, 

Dc(^)  =  -  0^(1  -  e  (4.29) 

results  in, 

e(>l»)  =  g^D^a^  +  gjg2aj^Dj(l  -  e  (4.30) 

where  H)  is  the  reduced  time  defined  in  Eq.  (3.2). 

Comparing  Eq.  (4.30)  with  (4.27), 

e  =  g_D  a_  (4.31) 

0  0  0  0 

ei(i')  =  giOj(l  -  e  (4.32) 

For  a  given  applied  stress  o^,  stress  developed  in  the  nonlinear  spring 
with  compliance  D^g^  is  and  the  corresponding  strain  is  e^.  For  the 
spring  with  the  nonlinear  compliance  D^gi,  the  stress  is  given  by. 


where  the  superscript  's'  denotes  quantities  related  to  the  springs. 
From  Fig.  1  it  is  evident  that  e|  and  e^^  are  equivalent.  Hence, 


al  =  g2(l  -  e  ^  )a^ 


(4.34) 


The  total  energy,  W^,  stored  in  the  two  springs  over  time  t  is  (see  Hiel 


et  al.  [551), 


0  t 

=  f  o^de  +  /  a^e^dt 
0  °  0  ^  ^ 


(4.35) 


Using  results  from  Eqs.  (4.31),  (4.33),  and  (4.34),  Eq.  (4.35)  becomes, 

'  1  Vo»o  "  ?  (4.36) 

For  a  viscoelastic  material  represented  by  multiple  Kelvin  elements  in 


series,  Eq.  (4.36)  takes  the  form. 


According  to  the  Reiner-Weisenberg  hypothesis,  failure  occurs  when  the 
stored  energy  reaches  the  resilience  of  the  material.  Denoting  the 
resilience  as  R,  the  expression  for  the  time  dependent  failure  stress 
obtained  from  Eq.  (4.37)  for  uniaxial  stress  state  is. 


(?  Vo  *  1  I  l“r(‘  - 


“  ^  -  'V  p 

!  "  )^ll 


(4.38) 
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4.4  Delayed  Failure:  Multiaxial  Formulation 

If  and  principal  stresses  at  any  point  in  a 

viscoelastic  material,  then  by  definition,  the  shear  stresses  are  zero 
on  the  principal  planes.  In  order  to  simplify  the  derivation,  let  it  be 
assumed  that  the  viscoelastic  material  is  represented  by  means  of  a 
single  Kelvin  element  (see  Fig.  1)  in  each  principal  direction.  The 
applied  multiaxial  creep  stresses  in  the  material  principal  directions 
are  given  by. 


=  a^H(t) 

(4.39) 

22  =  <^2”(t) 

(4.40) 

33  ^ 

(4.41) 

Substituting  Eqs.  (4.39),  (4.40)  and  (4.41)  in  Eqs.  (3.32),  (3.33)  and 
(3.35)  result  in  the  following  expressions  for  the  corresponding 
viscoelastic  strains, 

eil(t)  =  *  (1  -  *  (1  -  ^)03l  .  0,1(1  - 

0  0 

-X  1|)  J,  -n  -X  I) 

M(1  -  e  )  t  pi  (1  -  e  ))gjOj  *  (d  -  ?  ^  ) 

Ji 

+  ^  (1  -  ®  )}g2a3]  (4.42) 

From  Eq.  (4.42)  it  is  evident  that  the  effective  stress  developed  in  the 
spring  with  compliance  Dg  acting  in  principal  direction  1  is  given  by, 

J  J 

'’ol  ^  '^1  +  (1  -  o^)°3  (4.43) 

0  0 


■  ■.*>:*»*  ^ 4  j*4_a '  1^* 


Similarly,  the  effective  stress  developed  in  the  spring  with  compliance 


Oj^  and  acting  in  principal  direction  1  is. 


-X  i|)  -X  0)  J- 

(t)  =  (1  -  e  jgjOj  +  {(1  -  e  +  ji  (1  -  e  ^  ) 


'92»2 


-X  ii  J, 

+  {(1  -  e  *  ni  (1  -  e  ))g,0, 


(4.44) 


On  the  left  hand  side  of  Eqs.  (4.43)  and  (4.44),  the  superscript  's' 


denotes  the  effective  stress  within  the  spring,  the  first  subscript 


indicates  the  spring  number,  and  the  second  subscript  determines  the 


principal  direction  in  which  the  effective  stress  acts. 


The  total  energy,  W^,  stored  in  the  two  springs  in  material 
direction  1  over  time  t,  can  now  be  obtained  by  using  Eq.  (4.35). 


Hence, 


yS  _  1  n  /  S  \2  ^  1  n  /  ^ 


2 


(4.45) 


Using  a  procedure  similar  to  the  one  just  described,  it  can  be  shown 


that  for  an  isotropic  material  the  total  stored  energies  and  may 


be  expressed  in  a  form  similar  to  Eq.  (4.45).  Therefore,  the  total 


energy,  W^,  stored  in  the  springs  in  direction  j,  over  time  t  is  given 


“j  =  i  *  1  “I'-’u)' 


(4.46) 


If  the  viscoelastic  material  is  represented  by  n  Kelvin  units  in  series 


in  each  material  principal  direction  respectively,  then. 


is 


where. 


“j  -  i  “o(»oJ>^  ^  I  l,  ^(»rj)' 


s  /,  '^0./  V  ‘^0 

‘^Oj  "  ^  ‘^2  *  *C  ^i'lj 


(4.46) 


(4.47) 


.  J  -n  il» 

4i  =  92^^^  -  ®  )  -  r  ®  ''  )it°i 


+  02  +  0^\ 


92  Q  (1  -  e  i.j  =  1.2.3 


(4.48) 


Note  that  in  Eqs.  (4.47)  and  (4.48),  repeated  indices  imply  summation, 
and  6^j  is  the  Kronecker  delta  operator.  Also,  the  Prony  series  for  the 
creep  and  shear  compliance  are  required  to  have  the  same  number  of 
terms. 

Equations  (4.46),  (4.47)  and  (4.48)  define  the  energy  stored  in  the 
principal  direction  in  an  isotropic  viscoelastic  material. 

Therefore,  according  to  the  Reiner-Wei senberg  failure  theory,  the 
criterion  for  creep  rupture  in  the  jth  principal  direction  is  given  as. 


>  R 
J 

where  R  is  the  resilience  of  the  isotropic  material. 
For  a  material  with  a  constant  Poisson's  ratio, 

J(ii;)  *  (I  +  v)D(ij;) 


(4.49) 


(4.50) 


For  such  cases,  Eq.  (4.46)  is  still  valid,  but  Eqs.  (4.47)  and  (4,48) 
simplify  to. 


'  •  •  .  •  *  S.-  ’ 


CHAPTER  5 
NUMERICAL  RESULTS 
5.2  Preliminary  Comments 

In  this  section  results  of  a  number  of  linear  elastic,  linear 
viscoelastic  and  nonlinear  viscoelastic  analyses  are  discussed  in  light 
of  available  experimental  or  analytical  results.  All  results  are 
obtained  using  NOVA  on  an  IBM  3090  computer  in  double  precision 
arithmetic.  The  first  problem  deals  with  linear  elastic  (both  adhesive 
and  adherend)  analysis  to  show  the  effect  of  boundary  conditions  and 
mesh  on  the  stress  distributions.  Next,  results  of  geometric  nonlinear 
analysis  are  presented  and  compared  with  those  obtained  with  VISTA. 

Then  linear  and  nonlinear  viscoelastic  analysis  results  are  presented, 
first,  to  validate  the  finite  element  procedure  described  in  the 
preceeding  chapters  and,  second,  to  obtain  new  results  for  certain 
adhesive  joints. 


5.2  Linear  Elastic  Analysis :  Effects  of  Boundary 
Conditions  and  Mesh 

To  investigate  the  influence  of  boundary  cobnditions  on  the  elastic 
stress  distribution  in  a  single  lap  joint,  the  three  different  boundary 
conditions  shown  in  Fig.  2  were  used  in  the  linear  elastic  analysis. 
During  the  present  study  it  was  also  observed  that  the  type  of  finite- 
element  mesh  (i.e.  uniform  or  nonuniform)  has  also  an  effect  on  the 
stress  distribution  in  the  bondline.  The  material  properties  used  are 
given  in  Table  1. 


V ‘iT - -  - '  ’ 


Figures  3  and  4  show  plots  of  the  peel  stress  obtained  by  uniform 
and  nonunifortn  meshes,  respectively,  along  the  center  of  the  bond 
line.  Boundary  conditions  of  Type  1  and  3  give  almost  the  same 
distribution  of  the  stress,  while  Type  2  differs  significantly  at  the 
edges  of  the  adhesive.  Stresses  obtained  with  Type  1  and  3  boundary 
conditions  exhibit  stress  distributions  that  are  almost  symmetric  about 
the  vertical  centerline  of  the  joint  (with  Type  3  being  the  most 
symmetric).  It  is  also  observed  that  the  distribution  is  not  quite 
smooth  when  a  uniform  mesh  is  used.  For  an  accurate  description  of  the 
stress  gradients  near  the  edges,  a  more  refined  mesh  than  that  used  at 
the  center  (i.e.,  nonuniform  mesh)  must  be  used.  This  observation  is 
supported  by  the  results  shown  in  Fig.  3. 

The  effect  of  boundary  conditions  (Type  1  to  3)  on  the  distribution 
of  the  peel  and  shear  stresses  along  the  upper  and  lower  bondlines 
(i.e.,  interface  bewteen  the  adhesive  and  adherend)  are  shown  in  Figs. 
5-8.  The  nonuniform  mesh  was  used  in  all  cases.  From  these  results  it 
is  clear  that  boundary  condition  of  Type  2  gives  significantly  different 
results  than  Type  1  or  3,  especially  near  the  edges. 

5,3  Geometric  Nonlinear  Analysis 

Next,  geometrically  nonlinear  analysis  of  a  bonded  lap  joint  was 
considered.  The  geometry  and  boundary  conditions  of  Type  2  shown  in 
Fig.  2  are  used.  The  material  constants  used  are  given  in  Table  2.  The 
present  nonlinear  elastic  analysis  results  are  compared  with  those 
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ootained  using  the  VISTA  program  (23|  In  Figs.  9  and  10.  The  results 
are  In  excellent  agreement. 

Next,  the  nonlinear  response  of  a  bonded  cantilever  plate  under 
distributed  transverse  loads  was  Investigated.  The  plate  geometry  and 
the  two  finite  element  meshes  used  are  shown  in  Fig.  11.  The  material 
properties  used  are  given  in  Table  3.  Both  the  adhesive  and  the 
adherends  were  assumed  to  be  linearly  elastic  and  isotropic. 

The  load  on  the  plate  was  increased  in  steps  until  a  fairly  large 
free-end  deflection  was  obtained.  For  the  present  analysis  the 
magnitude  of  the  deflection  was  over  50%  of  the  beam  length.  The 
resulting  load-deflection  curves  obtained  by  the  two  meshes  are  shown  in 
Fig.  12.  The  results  obtained  by  using  linear  analysis  is  also  plotted 
for  comparison  purposes.  Clearly,  the  nonlinear  analysis  predicts  a 
stiffer  response.  This  is  due  to  the  fact  that  the  large  transverse 
deflection  causes  a  bending-extension  coupling  which  results  in  an 
increase  in  the  flexural  stiffness  of  the  beam. 

Figure  13  shows  the  compressive  bending  stress  at  a  specified  point 
(near  the  fixed  end)  in  the  lower  adherend  plotted  against  applied  load 
for  the  two  different  meshes.  The  discrepancy  in  the  two  curves  is  due 
to  the  fact  that  the  axial  stress  values  for  one  curve  were  obtained  at 
an  x-location  slightly  different  from  the  other  curve.  The  flattening 
out  of  the  stress  curve  at  higher  loads  is  a  result  of  the  shortening  of 
the  moment  arm  due  to  extensive  bending  of  the  beam. 
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Figures  14-16  show  the  variation  of  the  flexural,  peel  and  shear 
stresses  in  the  lower  half  of  the  adhesive  layer  plotted  along  the  plate 
axis  for  two  different  meshes.  Adjacent  to  the  clamped  end,  there 
exists  a  narrow  region  where  both  the  flexural  and  peel  stresses  are 
tensile.  However,  as  one  moves  further  along  the  plate  length,  the 
flexural  stress  turns  compressive,  which  conforms  to  what  is  predicted 
by  the  elementary  plate  theory.  The  shear  stress  attains  its  maximum 
value  near  the  clamped  end  and  decreases  rapidly  as  one  moves  out 
towards  the  free  end.  All  three  stresses  vanish  at  the  free  end  of  the 
plate,  thus  satisfying  the  stress  free  boundary  condition. 

5.4  Linear  Viscoelastic  Analysis 

The  nonlinear  constitutive  law  due  to  Schapery  may  be  linearized  by 
assuming  that  the  nonlinearizing  parameters  gg,  g^,  and  ^2  *  value 

of  unity.  In  addition,  the  stress  dependent  part  of  the  exponent  in  the 
definition  of  the  shift  factor  is  set  to  zero.  Consequently,  the 
constitutive  law  reduces  to  the  superposition  integral  form  commonly 
used  to  describe  a  linear  viscoelastic  material. 

Two  test  cases  are  used  to  validate  the  linear  viscoelastic 
analysis  capability  implemented  in  the  present  finite  element  program 
named  NOVA.  In  the  first  case,  the  tensile  creep  strain  in  a  single 
eight  noded  quadrilateral  element  was  computed  for  both  the  plane  stress 
and  plane  strain  cases  using  the  program  NOVA.  The  results  were  then 
compared  to  the  analytical  solution  for  the  plane  strain  case  presented 
in  (57].  A  uniform  uniaxial  tensile  load  of  13.79  MPa  was  applied  on 
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the  test  specimen.  A  three-parameter  solid  model  was  used  to  represent 
the  tensile  compliance  of  the  adhesive.  The  following  time  dependent 
functions  were  used  in  (57)  to  represent  the  tensile  compliance  and  the 
Poisson's  ratio  for  FM-73M  at  72‘’C: 


*  2[l+v(t)T  (5.1) 

Approximating  the  Poisson's  ratio  with  the  elasticity  relation  gives. 


»(t)  ■  3^  *  - 


(5.2) 


where  G(t)  and  K(t)  are  the  shear  and  bulk  modulus  (mm/mm/MPA) 
respectively,  and  Jg,  are  the  shear  compliance  coefficients.  The 
analytical  solution  to  the  creep  problem  for  the  plane  strain  case  is 
given  in  (57)  as: 

e(t)  =  2.728  X  10"^  +  1.334  x  10'^  -  2.659  x  lO'"*  g-t/0.3921 

It  is  to  be  noted  that  for  the  three-parameter  solid  charac¬ 
terization  of  FM-73M  the  value  of  the  Poisson's  ratio  actually  increases 
with  time.  However,  in  the  present  analysis  the  Poisson's  ratio  is 
assumed  to  be  independent  of  time.  Hence  two  discrete  values  of  the 
Poisson's  ratio  are  used  to  match  the  exact  solution  for  few  initial 
time  steps  and  final  time  steps.  The  values  of  the  Poisson's  ratio 

chosen  for  this  purpose  are  v„  =  Lim  v(t)  =  0.32  and  v  =  Lim  v(t) 

t-*0  t*™ 

=  0.417.  Figure  17a  shows  the  creep  curve  for  \>  =  0.417  for  both  plane 
strain  and  plane  stress  finite-element  analyses.  As  expected,  the  plane 
strain  results  exhibit  close  agreement  with  the  exact  solution  for  large 
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values  of  time,  followed  by  progressive  deterioration  of  predicted  value 
as  one  moves  towards  smaller  values  of  time.  The  finite  element  results 
for  the  plane  stress  case  points  to  the  fact  that  the  strains  are  higher 
for  plane  stress  than  for  plane  strain. 

Figure  17b  shows  the  creep  curve  corresponding  to  v  =  0.32  for  the 
plane  strain  case.  In  this  case  the  finite  element  predictions  are 
accurate  only  for  first  few  time  steps  and  deviates  more  and  more  from 
the  analytical  solution  as  time  increases.  This  is  not  surprising  since 
the  choice  of  Poisson's  ratio  for  this  case  makes  the  comparison 
meaningful  only  when  t  is  small. 

The  above  results  indicate  that  the  program  NOVA  provides 
reasonably  accurate  results  in  regions  where  the  input  parameters  are 
accurate,  and  that  the  variation  of  Poisson's  ratio  during  the  period  of 
analysis  may  cause  significant  deviations  from  the  actual  solution. 

Next,  the  Model  Joint  analysis  problem  presented  in  [57]  was  used 
as  the  second  validation  example.  In  this  case,  a  linear  viscoelastic 
finite  element  analysis  was  carried  out  on  a  model  joint  under  a 
constant  applied  load  of  4448  N  giving  an  average  adhesive  shear  stress 
of  13.79  MPa.  The  specimen  geometry,  discretization  and  boundary 
conditions  are  shown  in  Fig.  18.  The  thickness  of  the  adhesive  layer  is 
taken  to  be  0.254  mm.  A  nine  parameter  solid  model  was  used  to 
represent  the  tensile  creep  compliance  of  FM-73  at  72“C  and  is  given  by: 


D(t)  =  0.5988  X  10“^  +  1.637  x  10"®  (1  - 

+  0.6031  X  lO’"^  (1  - 
+  0.9108  X  10“^  (1  - 
+  2.6177  X  10"^  (1  - 

The  adhesive  Poisson's  ratio  is  assumed  to  have  a  value  of  0.417  and 
remains  constant  with  time.  The  material  properties  for  the  aluminum 
adherends  are  presented  in  Table  3. 

Figures  19  and  20  contain  plots  of  the  bond  normal  and  shear 
stresses,  respectively  for  t  =  50  secs,  and  t  =  60  min.  of  loading. 

These  stresses  represent  the  value  at  1/16  the  thickness  from  the  upper 
adhesive  adherend  interface.  The  sharp  peak  at  the  left  hand  edge  is 
due  to  the  singularity  caused  by  the  presence  of  a  re-entrant  corner  in 
the  vicinity  of  the  edge.  These  results  are  in  good  agreement  with  the 
results  presented  in  [571  which  uses  the  linear  viscoelastic  finite 
element  code,  MARC. 

5.5  Axisymmetric  Analysis  of  a  Linearly  Viscoelastic  Rod 
The  axial  displacement  of  one  end  of  a  linearly  viscoelastic  rod, 
subjected  to  a  spatially  uniform  end  traction  that  varies  sinusoidally 
with  time,  was  obtained  by  using  the  program  NOVA.  The  shift  factor  for 
the  material  is  defined  by  the  WLF  equation  and  the  temperature  is  held 
at  a  constant  value.  The  specimen  geometry  and  material  properties  are 


presented  in  Table  4.  The  exact  solution  to  this  problem  has  been 
presented  in  [231  and  was  used  to  validate  the  finite  element 
predictions.  As  can  be  seen  from  Fig.  21,  the  finite  element  results 
are  in  excellent  agreement  with  the  closed  form  solution  over  one  cycle 
of  loading  and  unloading. 

5.6  Nonlinear  Viscoelastic  Analysis  of  Adhesive  Coupons 

In  order  to  validate  the  nonlinear  viscoelastic  model,  three 
uniaxial  test  cases  are  analy2ed.  The  results  are  compared  with  the 
laboratory  tests  conducted  on  similar  specimens  by  Peretz  and  Weitsman 
[261.  The  material  properties  used  in  the  verification  analysis  are 
those  reported  in  [221.  The  creep  data,  together  with  other  relevant 
material  properties,  are  given  in  Table  5.  A  constant  value  for  the 
Poisson  ratio  is  assumed  for  the  adhesive.  The  results  from  a  linear 
viscoelastic  analysis  are  also  presented  for  comparison. 

In  the  first  verification  test,  a  uniaxial  stress  of  10  MPa  is 
applied  to  the  adhesive  coupon  for  1200  secs.,  followed  by  a  step 
increase  to  26.6  MPa  for  a  further  1200  secs.  The  temperature  of  the 
specimen  is  held  constant  at  50“C  and  is  assumed  to  be  uniform 
everywhere.  The  finite  element  predictions  for  this  test  are  plotted 
together  with  the  experimental  data  in  Fig.  22.  The  predictions  are  in 
good  agreement  with  the  experimental  results  of  Peretz  and  Weitsman 
[261. 

The  second  test  involves  creep  predictions  under  simultaneously 
varying  stress  and  temperature,  both  increasing  linearly  with  time.  The 
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temperature  1s  again  assumed  to  be  uniform  throughout  the  test 
specimen.  The  finite  element  predictions  (linear  and  nonlinear)  and 
[  experimental  data  are  compared  in  Fig.  23.  There  is  a  good  agreement 

j  between  the  two  sets  of  results 

J  The  third  test  involves  creep  under  a  constant  stress  of  10  MPa 

f 

<  with  a  linearly  varying  temperature  as  a  function  of  time.  Figure  24 

shows  the  strain  vs.  time  curves  obtained  in  the  experiments  and  finite 
element  analysis.  Satisfactory  agreement  between  the  experimental 
results  and  the  analysis  is  observed. 

A  further  set  of  tests  were  conducted  in  order  to  evaluate  the 
accuracy  of  the  finite  element  code  for  the  case  where  creep  is  followed 
by  creep  recovery.  A  qualitative  depiction  of  the  loading  and  the 
resulting  creep  strain  is  given  in  Fig.  25.  Rochefort  and  Brinson  [61] 
presented  experimental  data  and  analytical  predictions  on  the  creep  and 
creep  recovery  characteristics  of  FM-73  adhesive  at  constant 
temperature.  The  Schapery  parameters  necessary  to  characterize  the 
viscoelastic  response  of  FM-73  at  a  fixed  temperature  of  SO^C  are 


obtained  by  applying  a  least  squares  curve  fit  to  the  data  presented  in 
[611.  The  resulting  analytical  expressions  for  the  creep  compliance 
function  D(ii;),  the  shift  function  a^,  and  the  nonlinear  parameters  g^, 
gj^  and  g2  are  presented  in  Table  6.  From  the  point  of  view  of 
programming  convenience  it  is  more  suitable  to  work  with  an  exponential 
series  than  a  power  law.  Hence  the  power  law  creep  compliance  function 
was  converted  to  an  equivalent  five  term  exponential  series  of  the  form 


given  by  Eq.  (3.5).  The  five  constant  coefficents  for  this  series  were 
obtained  by  means  of  fitting  a  curve  to  the  aforementioned  power  law 
function  and  then  minimizing  the  error  in  a  least-squares  sense.  The 
exponential  series  form  of  the  compliance  function  is  presented  in  Table 
7  and  it  is  plotted  against  the  power  law  curve  in  Fig.  26  for 
comparison. 

Figure  27  shows  the  geometry  of  the  tensile  dogbone  specimen  used 
to  carry  out  the  creep  and  creep  recovery  tests.  This  geometry  is 
identical  to  the  one  used  by  Rochefort  and  Brinson  [611.  Due  to  the 
symmetry  of  specimen  geometry  and  applied  load,  only  the  upper  right 
hand  quadrant  of  the  specimen  was  analyzed.  The  finite  element 
discretization  consists  of  two  elements  along  the  length  of  the  specimen 
and  one  element  in  the  width  direction.  Eight-node  quadrilateral  plane 
stress  elements  are  used  for  this  analysis.  A  constant  tensile  load  is 
applied  on  the  specimen  for  the  first  30  min.  followed  by  creep  recovery 
over  an  equal  length  of  time.  The  procedure  is  repeated  for  three 
different  stress  levels  at  a  fixed  temperature  of  30°C. 

The  stress  input  for  a  uniaxial  creep  and  creep  recovery  test  is 
given  by, 

a(t)  =  agH(t)  -  o^H(t  -  t^)  (5.3) 

where  H(t  -  t^)  is  the  unit  step  function,  and  tj^  is  the  time  at  which 
stress  is  removed. 

Substitution  of  Eq.  (5.3)  into  Eq.  (3.1)  coupled  with  a  power  law 
representation  for  the  compliance  yields. 


(5.4) 
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and 

v(t)  ■  ^  1(1  -  v>"  -  (V>"i  <5-5) 

for  the  creep  and  creep  recovery  strains  respectively.  In  the  above 
expression. 


is  a  nondlmensional  parameter,  and 
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represents  the  transient  component  of  creep  strain  just  prior  to 
unloading.  Hence,  Eqs.  (5.4)  to  (5.7)  provide  a  closed  form  solution  to 
Schapery's  nonlinear  single  integral  law  for  the  simple  load  history 
involving  creep  and  creep  recovery  given  by  Eq.  (5.3). 

Figures  28  a,  b,  and  c  show  the  results  of  the  finite  element 
analysis  plotted  along  with  the  curve  representing  the  closed  form 
analytical  solutions  for  applied  stress  levels  of  21,  17  and  14  MPa 
respectively.  The  finite  element  predictions  are  in  excellent  agreement 
with  the  closed  form  solutions  except  at  the  beginning  of  creep  and 
again  at  the  onset  of  creep  recovery.  This  discrepancy  is  clearly  due 
to  the  discrepancy  between  the  power  law  and  the  exponential  series 
representation  of  the  creep  compliance  function  AD(ii)),  as  shown  in  Fig. 
26.  The  presence  of  too  many  data  points  in  the  far  field  region  has 
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caused  the  least  square  curve  fit  to  give  less  weight  to  the  initial 
data  points  and  therefore  overlook  the  error  present  near  the  beginning 
of  the  time  axis.  The  complete  agreement  between  the  closed  form 
solution  and  the  finite  element  prediction  for  large  values  of  time 
corroborates  this  fact.  From  Fig.  28  it  is  also  evident  that  the  error 
in  the  predicted  value  of  strain  decreases  as  the  applied  stress  is 
reduced.  This  is  exactly  what  is  expected  since  the  stress  dependent 
nonlinear  parameters  g^  and  g2  act  as  scale  factors  on  the  transient 
component  of  the  creep  strain.  Thus,  a  reduction  in  the  applied  stress 
causes  the  values  of  gj^  and  g2  to  reduce,  which  results  in  a 
proportionate  reduction  in  the  error  magnitude. 

5.7  Linear  and  Nonlinear  Viscoelastic  Analysis  of  a  Model 
Joint 

The  loading,  boundary  conditions  and  specimen  geometry  used  in  this 
analysis  is  the  same  as  the  one  used  in  the  earlier  model  joint  (see 
Fig.  18).  In  addition,  the  same  nine  parameter  solid  model  was  used  in 
this  analysis.  A  linear  viscoelastic  finite  element  analysis  was 
carried  out  over  a  period  of  one  hour  at  a  constant  applied  load  of 
3336  M.  The  results  for  the  linear  analysis  are  shown  in  Figs.  29-30. 
The  sharp  peak  at  the  left  hand  edge  is  due  to  the  singularity  caused  by 
the  presence  of  a  re-entrant  corner  and  dissimilar  materials.  All 
stress  plots  show  the  same  basic  trend  in  that  the  stresses  are 
attempting  to  redistribute  themselves  to  achieve  a  more  uniform 
distribution. 
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For  the  nonlinear  viscoelastic  analysis  of  the  model  joint,  the 
same  specimen  geometry  and  material  properties  were  employed.  However, 
the  nonlinearizing  parameters  and  the  shift  function  were  no  longer  held 
constant,  but  were  allowed  to  change  with  the  current  stress  state 
within  the  adhesive  layer.  The  results  from  this  analysis  are  presented 
In  Figs.  31  and  32.  It  Is  Immediately  apparent  that  the  effect  of  the 
nonlinearity  causes  a  'softening'  of  the  adhesive,  leading  to  a  response 
that  Is  less  stiff  compared  to  the  linear  case.  Hence,  even  though  the 
applied  load  Is  the  same,  the  shearing  strain  for  the  nonlinear  case  is 
significantly  larger  as  compared  to  the  linear  case  (Figs.  30  and  32). 
Moreover,  the  Increment  In  creep  strain  for  the  nonlinear  case  Is  0.0058 
as  compared  to  0.0041  for  the  linear  case  over  the  same  period  of 
time.  This  Is  exactly  what  is  expected  since  the  nonlinear  model  takes 
Into  account  the  acceleration  of  creep  caused  by  the  stresses  within  the 
adhesive. 

The  effect  of  the  nonlinearity  on  the  stress  curves  (Figs.  29  and 
31)  Is  to  create  a  more  uniform  stress  distribution  by  reducing  the 
stress  peaks  near  the  edges  while  increasing  the  stresses  at  the  mid¬ 
section  of  the  overlap.  The  significant  reduction  of  the  stress  peaks 
effected  by  the  nonlinear  model  Is  very  important  from  a  design  point  of 
view  since  the  reduction  of  stress  levels  at  the  critically  stressed 
regions  results  In  an  Improved  joint  efficiency. 
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5.8  Elastic  Analysis  of  a  Composite  Single  Lap  Joint 

Renton  and  Vinson  [37]  used  a  closed  form  elastic  solution  to 
conduct  a  parametric  study  of  the  effect  of  adherend  properties  on  the 
peak  stresses  within  the  adhesive  in  a  composite  single  lap  joint.  A 
similar  parametric  study  was  carried  out  using  the  finite  element 
program  NOVA.  The  geometry,  finite  element  discretization  and  boundary 
conditions  for  the  composite  lap  joint  are  shown  in  Fig.  33.  For  the 
sake  of  simplicity,  only  identical  adherends  are  considered.  Each 
adherend  is  made  up  of  seven  laminas  of  equal  thickness.  The 
orthotropic  material  properties  for  a  lamina  are  given  in  Table  8.  In 
order  to  maintain  material  symmetry  about  the  laminate  mid-plane  and 
thus  eliminate  bending-stretching  coupling,  a  eo/Oo/-0o/Oo/-9o/Oo/9o  ply 
orientation  was  selected  for  the  analysis.  Note  that  this  type  of  ply 
orientation  places  the  9°  ply  iimnediately  adjacent  to  the  adhesive 
layer.  The  adhesive  used  is  FM-73  and  its  isotropic  linear  elastic 
properties  are  listed  in  Table  9.  The  adhesive  layer  is  modeled  using 
sixteen  eight-noded  quadrilateral  elements  along  its  length  and  two 
elements  through  its  thickness.  A  series  of  elastic  finite  element 
analyses  is  performed  to  study  the  effect  of  ply  orientation,  lamina 
primary  modulus  (Qii)f  and  geometric  nonlinearity  on  the  peak  stresses 
in  the  adhesive. 

In  order  to  study  the  influence  of  ply  orientation  on  the  adhesive 
stress  distribution,  stress  analyses  were  performed  for  9  =0°,  9  =  15’, 
9  =  45°,  and  9  =  90’  respectively.  The  results  are  shown  in  Figs.  34 


I 


UH.i' 


•  • 


1 1 1  I  I  t  Mri  I  f  frnm  i  “i  rr 


mmmmmsm&mmK 


;s 

i 


78 


and  35.  The  plots  show  the  variation  of  stresses  along  the  upper 
bondline  of  the  overlap.  The  parameter  x/c  is  the  normalized  distance 
from  the  bond  centerline  such  that  the  value  x/c  =  -1  corresponds  to  the 
left-hand  free  edge  of  the  bond  overlap.  It  is  evident  from  these 
figures  that  an  increase  in  the  ply  orientation  angle  9,  causes  the  peak 
stresses  to  increase  near  the  free  edge  of  the  bond  overlap.  The 
adherend  with  a  O^/QO”  ply  orientation  (cross-ply)  shows  a  28%  increase 
in  peel  stress  and  a  17*  increase  in  shear  stress  over  the  corresponding 
values  for  a  0“  (unidirectional)  ply  orientation.  This  is  not 
surprising  since  a  cross-ply  adherend  has  a  lower  bending  stiffness 
which  results  in  a  larger  lateral  deflection  causing  higher  stress 
concentrations  at  che  overlap  ends. 

The  influence  of  the  lamina  primary  modulus  (Q^)  on  adhesive  peel 
and  shear  stresses  can  be  seen  in  Figs.  36  and  37  respectively.  A  0° 
(unidirectional)  adherend  ply  orientation  is  used  for  this  analysis. 

The  two  figures  show  a  significant  increase  in  the  peak  adhesive  stress 
as  the  value  of  decreases.  This  is  understandable  as  a  more 
flexible  adherend  would  undergo  larger  bending  and  hence  produce  higher 
stress  concentrations  at  the  overlap  ends. 

Harris  and  Adams  [651  conducted  large  displacement  finite  element 
analyses  on  a  single  lap  joint  with  aluminum  adherends  and  observed 
significant  reduction  in  peak  stresses  at  the  edge  of  the  adhesive  as 
compared  to  linear  results.  In  order  to  observe  the  effect  of  geometric 
nonlinearity  on  a  single  lap  joint  with  laminated  composite  adherends,  a 
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large  displacement  analysis  was  performed  using  the  program  NOVA.  Due 
to  its  greater  susceptibility  to  bending,  cross-ply  laminated  adherends 
were  used  for  this  analysis.  The  results  can  be  seen  in  Figs.  38  and 
39.  The  geometrically  nonlinear  analysis  results  in  a  30%  reduction  in 
the  peak  peel  stress  and  a  15%  reduction  in  the  peak  shear  stress.  The 
horizontal  shifting  of  the  nonlinear  curves  is  due  to  the  configuration 
coordinate  update  required  by  the  large  displacement  analysis. 

5.9  Nonlinear  Viscoelastic  Analysis  of  a  Composite  Single 
Lap  Joint 

A  nonlinear  viscoelastic  analysis  of  a  lap  joint  made  of  composite 
material  was  carried  out  over  a  time  period  of  forty  hours  using  NOVA. 
The  specimen  geometry  and  the  finite  element  discretization  are  the  same 
as  for  the  elastic  analysis  as  shown  in  Fig.  33.  However,  instead  of  a 
uniform  end  traction,  a  uniform  end  displacement  of  0.363  mm  is  applied 
to  the  end  of  the  joint  and  is  held  constant  with  time.  The  adherends 
are  made  of  symmetric  cross-ply  laminates  whose  properties  are  given  in 
Table  8,  while  the  adhesive  used  is  FM-73  and  its  creep  compliance  and 
Schapery  parameters  can  be  found  in  Table  5. 

Figures  40  and  41  show  the  variation  of  shear  stress  and  shear 
strain  respectively  across  the  entire  bond  length  over  a  period  of  40 
hours.  The  sharp  peak  on  the  left-hand  edge  is  due  to  the  presence  of  a 
re-entrant  corner  and  also  due  to  the  difference  in  material  properties. 
Figures  42  and  43  provide  a  close-up  view  of  the  shear  stress  and  strain 
gradients  at  the  free  edge.  As  might  be  expected,  the  shear  stress 
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undergoes  relaxation  which  results  in  a  36%  decrease  in  the  peak  value 
at  the  left  hand  edge.  The  stresses  have  been  normalized  with  respect 
to  an  average  shear  stress  value  of  4.5  N/mm^.  The  peak  shear  strain, 
however,  shows  an  increase  of  35%  over  the  same  period  of  time. 
Similarly,  Figs.  44  to  47  reveal  that  while  the  peak  values  of  the  peel 
and  axial  stresses  decrease  by  26%  and  32%  respectively,  the 
corresponding  strains  show  a  respective  increase  of  63%  and  6%.  The 
reason  that  the  strains  increase  with  time  even  though  the  joint  end 
deflection  remains  fixed,  is  because  the  adherends  are  modeled  as 
elastic  continuums.  As  the  stresses  in  the  adhesive  relax  with  time, 
the  elastic  adherends  deform  to  attain  a  new  equilibrium  configuration 
and  this  leads  to  an  altered  state  of  strain  within  the  adhesive. 

Hence,  it  is  very  important  that  the  elastic  nature  of  the  adherends  be 
taken  into  account  in  an  analysis.  Also,  the  significant  increase  in 
adhesive  strains  with  time  is  a  viscoelastic  phenomenon  and  therefore  it 
cannot  be  predicted  by  means  of  a  purely  elastic  analysis.  This  fact 
emphasizes  the  need  to  model  the  adhesive  layer  as  a  viscoelastic  medium 
in  order  to  be  able  to  predict  the  long-term  durability  of  a  bonded 
joint. 


5.20  Nonlinear  Fickean  Diffusion  in  Polystyrene 

In  order  to  validate  the  diffusion  model  implemented  in  NOVA  and 
discussed  in  Chapter  4,  results  from  a  nonlinear  diffusion  analysis 
presented  in  [66]  are  used.  The  test  problem  involves  unsteady  sorption 
of  a  penetrant  in  a  semi-infinite  medium  for  a  diffusion  coefficient 
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that  is  an  exponential  function  of  penetrant  concentration,  that  is, 

0  =  Dg  exp  (kC/Cg).  Finite  element  predictions  were  obtained  for  k  = 
0.614  and  k  =  3.912  and  the  results  were  compared  with  the  published 
results  represented  by  the  solid  lines  in  Fig.  48.  Excellent  agreement 
is  observed  for  the  two  values  of  the  coefficient  k. 

Levita  and  Smith  (67]  conducted  experiments  to  study  gas  transport 
in  polystyrene  and  found  that  the  diffusion  coefficients  for  gases 
decreased  with  time  when  the  polystyrene  film  was  subject  to  a  constant 
uniaxial  strain.  This  effect  was  attributed  to  the  continuous  free 
volume  recovery  (densif ication)  in  the  polystyrene  specimen  at  constant 
strain.  The  study  also  indicated  that  larger  free  volume  elements 
decrease  in  size  faster  than  the  smaller  ones  as  volume  recovery 
progresses.  Using  the  results  published  in  (67]  as  a  guideline,  NOVA 
was  used  to  study  the  time  dependence  of  the  diffusion  coefficient  for 
carbon-dioxide  gas  in  a  polystyrene  film  at  constant  strain.  For  this 
case,  the  temprature  and  moisture  concentration  effects  presented  in  Eq. 
(4.9)  were  neglected,  resulting  in  a  diffusion  coefficient  that  is 
solely  a  function  of  the  transient  component  of  the  dilatational  strain 
which,  in  turn,  is  a  measure  of  the  change  in  the  free  volume.  Figure 
49  shows  the  variation  of  the  diffusion  coefficient  with  time  for  three 
different  strain  levels.  The  material  properties  for  polystyrene  which 
were  obtained  from  [681  are  given  in  Table  10.  From  Fig.  49  it  is 
evident  that,  independent  of  the  strain  level,  the  diffusion  coefficient 
reaches  a  peak  value  t  =  1  at  hour  and  then  slowly  decays  to  the 
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reference  value,  Dg.  This  behavior  can  be  attributed  to  an  initial 
increase  in  free  volume  due  to  the  application  of  the  uniaxial  strain, 
followed  by  a  continuous  recovery  in  free  volume  (densification)  at  a 
constant  strain  as  the  polystyrene  film  undergoes  relaxation.  A  larger 
applied  strain  produces  larger  initial  dilatation,  and  this  results  in  a 
higher  peak  in  the  diffusion  coefficient.  Figure  49  also  reveals  that 
the  time  rate  of  free  volume  recovery,  and  hence  the  time  rate  of 
decrease  in  the  diffusion  coefficient,  is  proportional  to  the  applied 
strain  level. 

The  influence  of  penetrant  molecule  size  on  the  diffusion 
coefficient  for  gases  in  polystyrene  was  studied  by  varying  the 
magnitude  of  the  material  parameter  B  in  Eq.  (4.9).  The  temperature  and 
strain  were  held  constant  at  SO^C  and  1.8%  respectively.  The  prediction 
obtained  from  NOVA  are  shown  in  Fig.  50  for  two  values  of  B.  The  faster 
rate  of  decrease  in  the  diffusion  coefficient  for  a  higher  value  of  8 
implies  that  the  larger  free-volume  elements  decrease  in  size  faster 
than  the  smaller  ones  as  volume  recovery  progresses.  The  NOVA 
predictions  are  qualitatively  in  good  agreement  with  the  resu.ts 
presented  in  [671. 

When  a  polymeric  material  is  in  the  rubbery  state,  equilibrium  is 
reached  very  rapidly  in  response  to  variations  in  temperature,  stress 
and  penetrant  conentration.  By  contrast,  a  material  in  the  glassy  state 
is  not  in  thermodynamic  equilibrium  and  the  response  of  the  free  volume 
to  changes  in  external  conditions  is  delayed.  This  metastable  state 
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causes  the  free  volume  to  slowly  collapse  with  time  until  equilibrium  is 
reached.  This  phenomenon  is  known  as  physical  aging  and  causes 
relaxation  processes  to  take  place  over  a  longer  time.  Struik  [69] 
proposed  that  for  a  material  in  the  glassy  state,  effective  time  x  is 


related  to  actual  time  t  by. 


t  t^  u 

^  =  r  (r-fr^  • 

0  ^e  ^ 


(5.8) 


where  tg  is  the  aging  time  at  the  start  of  service  life  or  testing 
and  u  is  a  constant  such  that  0  <  u  <  1.  For  such  a  material,  the 
definition  of  reduced  time  given  by  Eq.  (3.2)  is  no  longer  valid  and 


should  be  modified  to. 


=  f  (a!T)'^ds 


(5.9) 


where  a  %  is  the  shift  factor. 


The  effect  of  physical  aging  on  the  diffusion  coefficient  for 
carbon-dioxide  gas  in  polystyrene  was  studied  by  implementing  Eqs.  (5.8) 
and  (5.9)  in  NOVA.  The  values  of  temperature,  strain  and  tg  were  set  at 
50'’C,  1.8%  and  24  hours  respectively.  Figure  51  shows  that  an  increased 
physical  aging  denoted  by  a  higher  value  of  the  parameter  u,  causes  the 
diffusion  coefficient  to  decay  slower  than  the  one  for  which  u  is 
lower.  This  behavior  is  expected  since  increased  physical  aging  causes 
the  free  volume  recovery  to  take  place  over  a  longer  period  of  time. 

Note  that  when  there  is  no  physical  aging,  u  and  tg  are  equal  to  zero 
and  X  is  identically  equal  to  t. 


A 


I  a- ^  I  i  igj  «*il**»« 


5.11  Linear  Elastic  Analysis  of  a  Butt  Joint 

Alvazzadeh  et  a1.  (311  used  special  linear  elastic  interface 
elements  to  study  the  effect  of  adhesive  thickness  and  adhesive  Young's 
modulus  on  the  stresses  within  a  bonded  butt  joint.  A  similar 
parametric  study  was  carried  out  using  NOVA  where  both  the  adhesive  and 
adherend  were  assumed  to  be  linearly  elastic.  The  specimen  geometry  and 
loading  are  shown  in  Fig.  52.  Due  to  symmetry,  only  a  quarter  of  the 
butt  joint  was  modeled.  The  finite  element  discretization  is  shown  in 
Fig.  53,  together  with  the  boundary  conditions.  The  various  adhesive 
and  adherend  properties  used  in  the  parametric  study  are  given  in  Table 
11.  A  plane  stress  elastic  finite  element  analysis  was  performed  and 
the  normalized  shear  and  normal  stresses  plotted  along  the  interface 
close  to  the  free  edge.  Figures  54  and  55  show  the  influence  of  the 
ratio  b/e  (where  b  is  the  width  of  the  butt  joint  and  e  is  the  thickness 
of  the  adhesive  layer)  on  the  adhesive  shear  and  normal  stresses 
respectively.  It  is  observed  that  the  maximum  value  of  shear  stress  and 
the  minimum  value  of  normal  stress  are  nearly  equal  for  different  joint 
thicknesses.  The  influence  of  the  ratio  of  adhesive  to  adherend  Young's 
moduli  on  adhesive  stresses  are  shown  in  Figs.  56  and  57.  As  this  ratio 
increases,  the  maximum  shear  stress  and  the  maximum  normal  stress 
increase  in  value  for  b/e  =  60. 


5.12  Nonlinear  Viscoelastic  Analysis  of  a  Butt  Joint 
Including  Moisture  Diffusion 

The  effect  of  a  change  in  the  free  volume  of  a  polymer  on  its 
viscoelastic  response  was  discussed  by  Knauss  and  Emri  [27].  They  used 
the  unifying  concept  of  the  free  volume  by  considering  that  fractional 
free  volume  depends  on  three  variables:  temperature  T,  moisture 
concentration  c,  and  mechanically  induced  dilatation  9.  Lefebvre  et  al . 
[51]  extended  the  free  volume  concept  to  define  a  nonlinear  diffusion 
coefficient,  which  results  in  a  coupling  between  the  viscoelasticity  and 
the  diffusion  boundary  value  problems  (see  Section  4.1).  The  influence 
of  this  coupling  on  the  viscoelastic  response  and  moisture  diffusion 
within  the  adhesive  layer  of  a  butt  joint  was  investigated  by  using  the 
program  NOVA.  The  specimen  geometry  and  finite  element  discretization 
are  the  same  as  shown  in  Figs.  52  and  53,  respectively.  However, 
instead  of  a  uniform  end  traction,  a  uniform  axial  displacement  of  0.1 
mm  is  applied  at  the  end  of  the  joint  and  is  held  constant  with  time. 

The  adherends  are  made  of  aluminum  and  the  adhesive  used  is 
polystyrene.  The  various  material  properties  are  listed  in  Tables  10 
and  11.  The  selection  of  polystyrene  as  an  adhesive  was  prompted  by  the 
fact  that  it  is  one  of  the  few  polymeric  materials  that  have  their 
viscoelastic  properties  and  diffusion  parameters  adequately 
documented.  The  normalized  moisture  concentration  at  the  free  edge  of 
the  adhesive  layer  is  unity,  and  the  initial  concentration  throughout 
the  adhesive  layer  is  zero.  The  tests  are  conducted  at  the  reference 
temperature  of  SO^C. 
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Figure  58  shows  the  moisture  concentration  profiles  within  the 
adhesive  layer  at  three  different  times  when  there  Is  no  coupling.  In 
this  case  the  diffusion  coefficient  remains  constant  with  time,  that  Is, 
D  »  Dq.  Fig.  59  shows  the  moisture  concentration  profiles  for  the  case 
where  there  Is  viscoelastic  coupling  only,  that  Is,  when  the  diffusion 
coefficient  depends  on  the  transient  component  of  the  dllatatlonal 
strain.  Fig.  60  depicts  the  case  where  there  Is  full  coupling,  that  is, 
the  diffusion  coefficient  Is  a  function  of  the  dllatatlonal  strain  and 
the  moisture  concentration  at  any  given  point  In  the  adhesive. 
Conversely,  the  viscoelastic  shift  factor  Is  now  a  function  of  the 
dllatatlonal  strain  and  the  moisture  concentration  ^see  Eq.  4.12).  Fig. 
61  presents  the  results  for  each  of  these  three  cases  for  comparison  at 
time  t  »  8  hours.  From  these  figures  It  Is  evident  that  the  effect  of 
coupling  Is  to  accelerate  moisture  diffusion  in  the  adhesive.  The 
mechanically  induced  dilatation  together  with  the  swelling  due  to 
moisture  sorption  results  in  a  higher  free  volume  fraction  within  the 
adhesive  which,  according  to  Eq.  4.9,  causes  diffusion  to  proceed  faster 
over  the  same  period  of  time.  It  Is  to  be  noted  that  in  Fig.  61,  the 
curves  become  less  concave  as  the  coupling  increases,  which  Is  in  good 
agreement  with  the  results  published  In  [66]. 

Figures  62  to  65  show  the  variation  of  the  stresses  and  strains 
with  time  within  the  adhesive  layer  In  the  butt  joint  when  there  is  no 
coupling  due  to  moisture  induced  swelling.  Mathematically,  this  implies 
that  Y=  0  In  equations  4.9  and  4.12.  From  Figs.  62  and  63  it  is  evident 


that  the  stresses  do  not  relax  significantly  over  the  time  period  of  the 
analysis.  This  Is  because  the  elastic  adherend  acts  as  a  spring  causing 
the  adhesive  to  creep  even  though  the  Joint  end  displacement  remains 
fixed.  However,  there  Is  a  slight  relaxation  In  the  normal  stress  as 
one  moves  towards  the  center  of  the  bond.  The  large  increase  In  the 
strains,  as  shown  In  Figs.  64  and  65,  Is  due  to  the  creep  caused  by  the 
strain  recovery  in  the  elastic  adherend.  This  observation  Is  supported 
by  Fig.  66  which  shows  that  the  normal  strain  In  the  adherend 
Immediately  adjacent  to  the  interface  undergoes  significant  reduction 
with  time.  The  decrease  In  the  adherend  normal  stress,  as  shown  In  Fig. 
67,  reflects  the  concurrent  stress  relaxation  that  occurs  In  the 
adhesive  and  triggers  the  strain  recovery  In  the  adherend. 

Figs.  68  to  71  show  the  effect  of  moisture  Induced  swelling  on  the 
viscoelastic  stresses  and  strains  In  the  adhesive  layer. 

Mathematically,  this  means  y  has  a  nonzero  value  In  Eqs.  4.9  and  4.12. 
The  actual  value  of  y  selected  for  this  study  is  0.001.  For  this  value 
of  y,  the  moisture  absorbed  causes  large  swelling  strains  within  the 
adhesive,  which  increase  in  magnitude  as  the  diffusion  progresses.  This 
moisture  Induced  swelling  strain  causes  a  reduction  in  the  mechanically 
Induced  normal  strain  and  hence  a  lower  value  for  the  normal  stress  In 
the  adhesive.  This  effect  can  be  observed  In  Fig.  68  where  progressive 
swelling  has  caused  a  25%  reduction  In  the  peak  normal  stress  over  a 
period  of  8  hours.  It  Is  Interesting  to  note  that  the  difference 
between  the  two  stress  curves  diminishes  as  one  moves  towards  the  center 


of  the  bond.  This  behavior  Is  expected  since  there  1s  very  little 
moisture  near  the  center  of  the  bond  and  so  the  stress  reduction  is 
primarily  due  to  viscoelastic  relaxation.  The  large  increase  in  the 
adhesive  strains,  as  seen  In  Figs.  70  and  71,  Is  due  to  the  adherend 
acting  as  a  elastic  spring. 

Fig.  72  shows  the  Influence  of  the  moisture  coefficient  y,  on  the 
normal  stress  In  the  adhesive  layer  after  eight  hours  of  sorption.  As 
can  be  seen,  the  swelling  Induced  for  y  *  0.001  results  In  a 
significantly  lower  normal  stress  near  the  free  edge  as  compared  to  the 
case  where  y  >  0.  Away  from  the  free  edge,  the  two  stress  curves  appear 
to  merge  as  one  moves  towards  the  Interior  of  the  bond.  This  Is  because 
the  low  moisture  concentrations  present  In  the  bond  Interior  Is 
Insufficient  to  cause  any  significant  reduction  In  the  normal  stress  due 
to  swelling. 

5.13  Delayed  Failure  of  a  Butt  Joint 

The  theory  presented  in  Secs.  4.3  and  4.4  was  applied  to  predict 
viscoelastic  creep  failure  within  the  adhesive  In  a  butt  joint.  The 
specimen  geometry  and  the  finite  element  discretization  are  the  same  as 
shown  In  Figs.  52  and  53,  respectively.  The  adherend  is  made  of 
aluminum  and  Its  material  properties  are  given  In  Table  11.  The 
adhesive  used  Is  FM-73  and  its  tensile  creep  compliance  is  listed  in 
Table  5.  The  failure  parameter  (R)  for  FM-73,  also  known  as  the 
resilience,  was  obtained  by  computing  the  area  under  the  stress-strain 
curve  presented  in  [70].  This  procedure  yielded  a  value  of  the 
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resilience  as  1.3  N.mm/mm^.  Note  that  the  area  under  the  visco-plastic 
yield  plateau  was  not  Included  In  computing  the  value  of  R.  According 
to  the  Relner-Welsenberg  theory,  failure  occurs  when  the  stored  energy 
per  unit  volunra  In  the  body  reaches  the  resilience  R,  for  the  material. 
Using  this  postulate  as  a  failure  criterion,  NOVA  was  utilized  to 
predict  failure  In  the  adhesive  layer  of  the  butt  joint  subject  to  a 
constant  uniaxial  tension.  The  Influence  of  applied  stress  level  on 
delayed  failure  was  studied  by  using  a  stress  level  of  69,  60,  and  54 
MPa  respectively.  In  all  three  cases,  failure  was  Initiated  In  the 
adhesive  element  located  right  at  the  free-edge  and  Immediately  adjacent 
to  the  Interface.  It  was  also  observed  that  the  direction  of  the  plane 
of  failure  was  always  Inclined  at  an  angle  of  18®,  counter-clockwise  to 
the  x-ax1s.  Since  the  direction  of  failure  coincides  with  the  direction 
of  principal  stress,  it  is  evident  that  a  multiaxlal  state  of  stress 
exists  near  the  free  edge,  even  though  the  applied  stress  is  uniaxial. 
This  observation  is  in  agreement  with  the  results  presented  in  Secs. 

5.11  and  5.12.  Fig.  73  shows  the  variation  of  normal  (or  creep)  strain 
with  time  at  30"C  for  the  element  In  which  failure  Is  first  initiated. 
The  right  hand  termination  point  on  the  curves  indicate  the  point  at 
which  failure  occured.  It  is  observed  that  for  an  applied  stress  level 
of  69  MPa,  the  time  to  failure  (tp)  Is  1.5  secs.  In  other  words,  for 
this  stress  level,  failure  occurs  almost  instantaneously.  For  an 
applied  stress  of  60  MPa,  tp  increases  to  400  secs.  Reducing  the 
applied  stress  to  54  MPa  results  in  a  time  to  failure  of  approximately 


10  hours.  These  results  are  qualitatively  in  good  agreement  with  the 
results  presented  by  Bruller  [45]  for  PMMA. 


ir{ 


From  the  above  observations  it  is  clearly  evident  that  for 
viscoelastic  polymers  like  FM-73,  the  time  to  failure  depends  strongly 
on  the  applied  stress  level.  Fig.  74  shows  the  evolution  of  stored 
energy  with  time  for  different  stress  levels.  For  very  high  applied 
stress  levels,  almost  all  the  strain  energy  is  conserved  as  stored 
energy  and  failure  occurs  almost  immediately.  For  intermediate  levels 
of  applied  stress,  viscoelastic  creep  causes  a  part  of  the  strain  energy 
to  be  dissipated.  As  a  result,  only  a  fraction  of  the  total  strain 
energy  is  conserved  as  stored  energy.  Consequently,  the  stored  energy 
builds  up  slowly,  analogous  to  a  "leaking  vessel",  resulting  in  delayed 
failure.  For  an  applied  stress  level  that  is  below  a  certain  threshold 
value  for  a  given  material,  the  dissipated  energy  may  constitute  a  large 
fraction  of  the  total  strain  energy.  In  that  case,  the  stored  energy 
would  increase  too  slowly  to  exceed  the  resilience  of  the  material  over 
any  realistic  length  of  time,  and  hence  there  would  be  no  failure  even 
if  the  applied  stress  acts  indefinitely. 
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CHAPTER  6 


SUMMARY  AND  CONCLUSIONS 

6.1  General  Summary 

A  nonlinear  viscoelastic  computational  model  Is  developed, 
validated  and  applied  to  the  stress  analysis  of  adhesively  bonded 
Joints.  The  large  displacements  and  rotations  experienced  by  the 
adherends  and  the  adhesive  are  taken  Into  account  by  Invoking  the 
updated  Lagranglan  description  of  motion.  The  adhesive  layer  Is  modeled 
using  Schapery's  nonlinear  single  Integral  constitutive  law  for  uniaxial 
and  multlaxial  stress  states.  The  effect  of  temperature  and  stress 
level  on  the  viscoelastic  response  Is  taken  into  account  by  a  nonlinear 
shift  factor  definition.  Optionally,  a  nonlinear  shift  factor 
definition  based  on  the  concept  of  free  volume  that  was  postulated  by 
Knauss  Is  also  available.  Penetrant  sorption  is  accounted  for  by  a 
nonlinear  Fickean  diffusion  model  In  which  the  diffusion  coefficient  is 
dependent  on  the  temperature,  penetrant  concentration,  and  the 
dllatatlonal  strain.  A  delayed  failure  criterion  based  on  the  Reiner- 
Weisenberg  failure  theory  has  also  been  implemented  in  the  finite 
element  code.  The  program  is  validated  by  comparing  the  present  results 
with  analytical  and  experimental  results  available  in  the  literature. 
Additional  results  for  a  bonded  cantilever  plate,  single  lap  joint, 
thick  adherend  specimen,  and  butt  joint  are  also  presented.  The  program 
capability  has  been  extended  to  account  for  laminated  composite 
adherends  and  adhesives  with  a  time  dependent  Poisson's  ratio.  In 
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general,  the  computer  program  developed  herein,  named  NOVA,  is  believed 
to  provide  accurate  predictions  over  a  wide  range  of  specimen 
geometries,  external  loads,  and  environmental  conditions. 


6 . 2  Concl usions 

The  results  presented  in  Ch.  5  underscore  the  importance  of 
modeling  the  adhesive  in  a  bonded  joint  as  a  viscoelastic  material. 

This  allows  the  analyst  to  predict  the  large  increments  in  adhesive 
strains  that  occur  with  time  and  cannot  be  predicted  by  a  purely  elastic 
analysis.  Furthermore,  other  events  (such  as  moisture  diffusion  and 
delayed  failure),  that  are  highly  relevant  for  bonded  joint  analysis, 
cannot  be  accurately  predicted  unless  viscoelasticity  is  taken  into 
account.  At  high  stress  levels,  nonlinear  viscoelastic  effects  can 
produce  creep  strains  that  are  significantly  larger  than  the  linear 
viscoelastic  predictions  and  such  effects,  therefore,  should  be 
accounted  for.  The  effect  of  change  in  Poisson's  ratio  with  time  in 
some  polymers  have  a  significant  bearing  on  the  final  response  and  must 
be  taken  into  account  in  order  to  obtain  accurate  results. 

The  results  in  Chapter  5  also  indicate  that  the  stress  boundary 
conditions  at  the  free  edges  of  the  adhesive  are  not  exactly 
satisfied.  This  deficiency  in  the  model  is  expected  because  a 
displacement  based  finite-element  formulation  satisfies  the  boundary 
conditions  only  in  a  global  sense.  Even  so,  the  shear  stress,  as 
presented  in  Chapter  5,  shows  a  tendency  to  drop  towards  zero  as  it 
approaches  the  free  edge.  Any  deviations  from  this  behavior  can  be 


Table  2.  Data  for  Geometric  Nonlinear  Analysis  of  a  Lap  Joint. 


'?■ 

I 

•r. 


Table  3.  Data  for  Geometric  Nonlinear  Analysis  of  a  Bonded  Cantilever 
Plate. 


Adherend  fAlumlnuml 


E  «  70  X  10^  MPa 


Adhesive 


E  =  2.8  X  lO-*  MPa 
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Table  5.  Material  Data  for  FM-73  Unscrlmmed  at  30°C. 


Elastic  Compliance,  0^: 

Poisson's  Ratio,  v; 

Coefficient  of  Thermal  Expansion,  a: 
Prony  Series  Coefficients: 


360  X  lO'^/MPa 
0.38 

6.6  X  10"^  m/m/'’IC 


0l  »  11.05xl0-®/MPa 
O2  *  12.27xl0-®/MPa 
D3  =  17.35xl0‘®/MPa 
D4  =  21.63xlO-®/MPa 
O5  =  31.13xlO-®/MPa 
Dq  =  41.78xlO-®/MPa 


Tj  »  10  secs. 

»  10^  secs. 

Tj  =  10^  secs. 

T-  =  10^  secs. 
4 

Tg  =  10®  secs. 
T,  =  10®  secs. 


Table  6.  Data  for  Creep  and  Recovery  of  FM-73  Adhesive. 

0(0)  *  Dq  '*■ 

Oq  »  227.573  X  10“®/MPa 

O^(tt))  =  C4>" 

c  »  31.763  X  10"®/MPa 
n  *  0.151 
a  »  1  -  3.536  X 

a 

»  1  +  2.247  X 

^0 

=  1  +  6.981  X  10"^a^*®® 
g2  =  1  +  3.098  x  10' 
where  o  Is  In  MPa. 


l*t  **t.i*t. 


Table  8.  Orthotropic  Material  Properltes  for  Composite  Adherend 


Qll  »  46.885x103  MPa 
Qi2  -  Qi3  »  4.137x103  MPa 
Q22  =  Q33  •  14.962x103  MPa 
Q23  *  ^32  *  2.068x10^  MPa 
Q44  *  Q55  •  ^66  *  3.447x10^  MPa 


Table  10.  Material  Properties  for  Polystyrene  at  50®C. 


Bulk  Compliance: 

Mg  »  1.2xlO“VMPa 
Mj  *  0.2896xl0-^/MPa 
M2  -  0.2246xl0"^/MPa 
M3  »  0.3721xl0-^/MPa 
M4  -  0.1354xl0-^/MPa 
Shear  Compliance: 

Jq  »  l.OxlO-3/MPa 
»  2.16/MPa 
J2  -  2.92/MPa 
J3  »  1.38/MPa 
J4  »  2.88/MPa 
J5  »  2.31/MPa 
Jg  =  3.59/MPa 
Jy  »  0.648/MPa 


"1 

^2 

"3 

"4 


'’I 

”2 

’’3 

’^4 

'^5 

'’6 

'’7 


Reference  free  volume  fg  »  0.033 
Diffusion  coefficient  Og  *  9x10"®  mm^/sec 


»  1.515x10^  sec. 

»  1.515x10^  sec. 

=  1.515x10^  sec. 

“  1.515x10®  sec. 

»  1.515x10®  sec. 

»  1.515xlolO  sec. 
»  1.515x10^^  sec. 
»  1.515x10^®  sec. 
»  1.515x10^^  sec. 
=  1.515x10^®  sec. 
=  1.515x10^®  sec. 


Table  11.  Properties  for  Elastic  Analysis  of  a  Butt  Joint. 


Materials 

E(MP») 

V 

Steel 

2,07x10^ 

0.29 

Aluminum 

0.7x105 

0.33 

Eponal 

5,8x103 

0.33 

Rigid  Epoxy 

2.2x103 

0.33 

Ilf 

1 

V 

i 

I 

I 

I 


.f* 

I 

J** 

V, 


Intensity  of  the  transverse  load. 


Figure  12.  Load  -  Deflection  Curres  for  a  Bonded  Cantilcrer  Plate. 


Oondline  for  a  Bonded  Cantilever  Plate. 


1,000 


Dondlinc  for  a  Bonded  Cantilever  Plate. 
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Figure  21.  Plot  of  Axial  Displacement  against  Time  for  a  Viscoelastic  Rod 
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Figure  22.  Variation  of  Creep  Strain  with  Time  for  a  FM-73  Coupon  Subject  to  Step 
Loads  at  Constant  Temperature. 


Stress  0^  (MPa) 


Percent  strain 


Figure  23.  Stress-Strain  Curre  for  a  FM-73  Coupon  Under  Linearly  Varying  Stress  and 
Temperature. 
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Figure  24.  Varladon  of  Creep  Strain  with  Time  for  a  FM-73  Coupon  Subject  to  a 
Linearly  Varying  Temperature  at  Constant  Stress. 
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Figure  34.  Inducnce  of  Ply  Orientation  on  AUhesive  Peel  Stress. 
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Figure  38,  Influence  of  Geometric  Nonlinearity  on  Ailhesive  PeeJ  Stress. 
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Figure  40.  Variation  of  Shear  Stress  with  Time  for  Entire  Overlap. 
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Figure  44.  Variation  of  Peel  Stress  with  Time  Near  the  Free  Edge. 
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Figure  46.  Variation  of  Axial  Stress  with  Time  Near  the  Free  Edge. 
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Figure  47.  Variation  of  Axial  Strain  with  lime  Near  the  Free  EJge. 
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Figure  57.  Influence  of  Modulus  Ratio  on  Adhesive  Nomtal  Stress. 
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Figure  64.  Variation  of  Normal  Strain  in  the  Adhesive  With  Time  For  Viscoelastic 
Coupling. 


KKTOXffVUtB 


0.8  -rr- 

0.92 


0.94  0.98 

DBTANCE,  Y/  Sayf^-m 


0.98 


1.00 


Figure  67.  Varution  of  Normal  Strain  in  the  Adhcrcnd  With  Time  For  Vlscociasiic 
Coupling. 
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Figure  68.  Variation  of  Normal  Stress  in  the  Adhesive  With  Time  For  Viscoelastic  and 
Moisture  Coupling. 
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Figure  69.  Variation  of  Shear  Stress  in  the  Adhesive  With  Time  For  Viscoelastic  and 
Moisture  Coupling. 
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Figure  70.  Variation  of  Normal  Strain  in  the  Adhesive  With  Time  For  Viscoelastic  and 
Moisture  Coupling. 
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Figure  72.  Influence  of  Moisture  Coefficienl  on  Adh«sj»r  ^  st-rv*  < 
Hours. 
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